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CHAPTER  I 


INTRODUCTION 


A.  Heterojunctions 

Semiconductor  heterojunctions  are  junctions  between  chemi- 
cally different  semiconductors,  such  as  between  GaAs  and  ALAs , or 
between  GaAs  and  Ge.  The  two  semiconductors  which  form  a 
heterojunction  may  have  significantly  different  properties,  so  the 
use  of  such  junctions  promises  to  increase  vastly  the  flexibility 
available  in  the  design  of  electronic  devices.  However,  the  tech- 
nology involved  in  fabricating  heterojunctions  is  complex,  and 
often  is  not  readily  transferable  from  one  material  to  another.  Thus 
there  is  a need  to  be  able  to  predict  the  properties  of  a hetero- 
junction  between  two  given  semiconductors  before  that  hetero- 
junction is  in  fact  fabricated.  In  this  work  we  describe  a theoretical 
model  which  successfully  gives  the  energy  band  lineups  of  known 
heterojunctions,  and  which  is  readily  applied  to  unknown  ones. 

Sharma  and  Purohit  have  published  an  exhaustive  review  of 
heterojunction  work  up  to  about  197  3 [Sharma,  1974],  so  only  a 
brief  historical  account  will  be  given  here.  The  first  proposal  to 
use  different  semiconductors  in  the  same  electronic  device  was 
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made  by  Shockley  in  the  patent  on  the  bipolar  junction  transistor 
[Shockley,  1951].  Real  interest  in  heterojunction  technology  did 
not  develop,  however,  until  Kroemer  showed  that  significant 
improvements  in  transistor  performance  could  be  achieved  by  using 
a wider-gap  semiconductor  for  the  emitter  than  that  of  the  base 
[Kroemer,  1957].  Anderson  fabricated  the  first  heterojunctions 
[Anderson,  1960]  and  presented  a theoretical  model  [Anderson, 

1962]  for  the  band  lineup  which  is  commonly  used  today  to  explain 
abrupt,  lattice-matched  heterojunctions.  Heterojunction  work  was 
further  stimulated  when  Kroemer  suggested  the  double-heterostructure 
laser  [Kroemer,  1963].  Currently  the  most  important  application  of 
heterojunction  technology  is  in  such  opto-electronic  devices,  and 
they  are  rapidly  approaching  commercial  production  [Kressel,  197  6]. 

The  Anderson  model  of  the  band  lineup  at  a heterojuncticn  is 
illustrated  in  Figure  1,  with  band  energy  plotted  as  a function  of 
position  through  the  junction.  Anderson  assumes  that  the  bulk 
semiconductor  band  structures  exist  up  to  the  junction  where  they 
discontinuously  change  from  one  to  the  other.  The  magnitude  of  the 
discontinuity  in  the  conduction  band  is  assumed  to  be  given  by  the 
difference  in  the  electron  affinities  between  the  two  semiconductors. 
The  electron  affinity  x is  defined  as  the  energy  required  to  move  an 
electron  from  the  bottom  of  the  conduction  band  to  a point  just  out- 
side the  semiconductor  (the  "vacuum  level").  Thus  the  basic 


assumption  in  the  Anderson  model  is  that  the  vacuum  levels  line  up 


ELECTRON  ENERGY 


Figure  1.  The  Anderson  model  for  the  energy  band  profile  of  an 


abrupt  hetero junction. 
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at  the  interface.  Once  the  heterojunction  is  formed,  the  Fermi 
levels  in  the  two  semiconductors  must  be  equalized.  This  involves 
a transfer  of  electrons  from  one  side  of  the  junction  to  the  other, 
creating  a space  charge  with  the  consequent  band  bending. 

The  Anderson  model  has  not  been  completely  successful  in 
explaining  the  transport  properties  of  experimental  heterojunctions, 
so  several  other  models  have  been  proposed  taking  into  account 
other  processes.  Dolega  assumed  that  there  are  a large  number  of 
localized  electron  states  at  the  interface  [Dolega , 1963].  Rediker, 
et  al.,[1964]  and  Riben  and  Feucht  [1966]  studied  various  tun- 
neling processes.  In  this  work  we  are  interested  in  the  electronic 
structure  of  an  ideal  heterojunction.  In  other  words,  we  wish  to 
know  the  energy  and  location  of  the  electronic  states  in  the  vicinity 
of  the  interface.  We  are  only  secondarily  interested  in  what  the 
electronic  structure  implies  about  he  transport  properties.  There- 
fore, we  assume  most  of  the  features  of  Anderson's  model:  discon- 
tinuous change  in  band  structure  at  the  interface  and  rigid  bands. 

Our  only  significant  departure  from  the  Anderson  model  has 
been  to  regard  the  conduction  band  discontinuity  as  a property  of 
the  interface,  to  be  calculated  with  a model  appropriate  to  hetero- 
junctions, rather  than  assuming  it  to  be  determined  by  the  difference 
in  electron  affinities,  as  Anderson  did.  Kroemer  has  discussed  the 
arguments  against  accepting  the  electron  affinity  rule  [Kroemer, 

1975].  Essentially  the  problem  is  that  the  electron  affinity  is 
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largely  dependent  on  the  surface  dipole  set  up  by  the  electron  dis- 
tribution at  a free  surface.  There  is  no  reason  to  believe  that  at  a 
heterojunction  these  dipoles  will  add  linearly,  since  the  electron 


distribution  is  almost  certain  to  resemble  more  nearly  the  bulk 
semiconductor  than  a superposition  of  free  surfaces.  Therefore, 
in  the  absence  of  information  on  the  electron  distribution,  the  elec- 
tron affinity  rule  is,  in  principle,  incapable  of  predicting  hetero- 
junction band  lineups. 

Our  approach  has  been  to  replace  the  electron  affinity  rule 
with  one  that  is  somewhat  similar,  but  involving  a different  quan- 
tity. Instead  of  expressing  the  energy  of  the  bands  in  terms  of  a 
vacuum  level  defined  at  a location  outside  of  the  crystal,  we  refer 
the  band  energies  to  a level  defined  by  the  electrostatic  potential 
within  the  crystal.  We  have  chosen  a reference  level  which  we 
have  reason  to  believe  is  continuous  across  most  heterojunctions 
of  interest.  Thus  the  problem  is  reduced  to  the  calculation  of  both 
the  band  structures  and  the  electrostatic  potential  in  terms  of  a 
common  energy  scale.  We  have  done  so  using  a self-consistent 
pseudopotential  technique  which  will  be  described  in  detail  later. 


B.  Pseudopotentials 

The  pseudopotential  method  for  band  structure  calculations  is 
the  foundation  upon  which  most  of  this  work  rests.  It  was  chosen 
because  of  the  relative  ease  with  which  one  can  obtain  results 
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applicable  to  real  semiconductors,  in  contrast  with  more  formal 
methods  like  the  orthogonalized  plane  wave  (OPW)  or  augmented 
plane  wave  (APW)  methods.  It  has  also  been  the  only  method  which 
has  so  far  proved  feasible  for  calculations  of  the  electronic  struc- 
ture of  surfaces,  so  it  would  be  the  obvious  choice  if  one  wished 
to  extend  this  work  to  detailed  calculations  of  the  interface. 

An  excellent  review  of  pseudopotential  work  is  available  in 
Heine  [1970]  and  Cohen  [1970]  so  again  only  a short  summary  will 
be  given  here.  The  motivation  for  the  pseudopotential  method  was 
the  difficulty  of  treating  the  extended  (valence  and  conduction  band) 
electron  states  with  the  same  set  of  basis  functions  as  those 
required  for  the  tightly  bound  core  states.  The  core  states  are  best 
described  in  terms  of  localized,  atomic-like  states,  while  the 
extended  states  are  best  described,  in  the  region  outside  the  core, 
by  plane  waves.  The  OPW  method  of  Herring  [1940]  tried  to  over- 
come these  difficulties  by  setting  up  the  calculation  on  a basis  of 
plane  waves  which  had  been  orthogonalized  with  respect  to  the  core 
states.  This  method  requires  rather  massive  computations,  so  its 
application  has  been  somewhat  limited. 

The  pseudopotential  concept  was  first  proposed  by  Phillips 
and  Kleinman  [1959].  They  noted  that  the  orthogonahzing  terms  in 
the  OPW  method  could  be  grouped  with  the  potential  in  the 
Schroedinger  equation,  giving  rise  to  an  effective  (=pseudo) 


potential 
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V , = V + Z (E-E  ) 1c)(cI  (I.  1) 

Ic) 

where  the  sum  is  over  all  core  states  Ic) . The  orthogonality  terms 
are  expected  to  have  the  effect  of  a repulsive  potential,  so  the 
pseudopotential  should  be  much  weaker  than  the  irue  potential. 
Cohen  and  Heine  [M.  H.  Cohen,  1961]  considered  the  Phillips- 
Kleinman  form  in  more  detail  and  showed  that  this  is  indeed  the 
case. 

Heine  and  Abarenkov  [Heine,  1964;  Abarenkov,  1965] 
developed  a model  pseudopotential  for  the  ion  core,  consisting  of 
the  Coulomb  potential  outside  of  a given  core  radius,  and  a con- 
stant potential  inside,  which  was  fitted  separately  for  each  angular 
momentum  value.  These  parameters  were  determined  by  fitting  the 
spectroscopically  observed  energy  levels  of  the  closed-shell-plus- 
one-electron  ion.  Animalu  and  Heine  [1965  and  1966]  determined 
the  parameters  for  a large  number  of  elements . These  pseudo- 
potentials  have  been  widely  used  as  initial  potentials  which  were 
then  varied  to  fit  empirical  data. 

Brust  [1964]  applied  the  empirical  pseudopotential  method 
(EPM) , in  which  the  pseudopotential  form  factors  are  treated  as 
freely  adjustable  parameters,  to  silicon  and  germanium.  He 
developed  the  calculational  techniques  which  are  still  commonly 
used.  This  involves  using  a basis  set  of  plane  waves  with  a cutoff 
at  a given  magnitude  of  wavevector  and  treating  other  plane  waves 
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by  a modification  of  the  Lbwdin  perturbation  theory  [Lbwdin,  1951]. 

Cohen  and  Bergstresser  [M  . L.  Cohen,  1966]  applied  the  EPM  to  the 
calculation  of  the  band  structures  of  fourteen  semiconductors, 
demonstrating  the  power  of  the  pseudopotential  method  in  dealing 
with  real  materials.  The  EPM  has  since  been  used  to  calculate  such 
diverse  properties  of  bulk  semiconductors  as  dielectric  functions 
[Richardson,  197 1;  Vinsome,  1971],  electronic  densities  of  states 
[Chelikowsky , 1974b],  and  the  temperature  dependence  of  X-ray 
reflection  intensities  [Chelikowsky,  1974a]. 

The  extension  of  pseudopotential  calculations  to  systems  more 
complicated  than  perfect  crystals  was  made  possible  by  the  self- 
consistent  pseudopotential  method,  first  suggested  byAppelbaum 
and  Hamann  [1973].  The  motivation  for  developing  this  technique 
was  the  study  of  the  electronic  structure  of  semiconductor  surfaces, 
and  in  this  application  it  has  been  quite  successful  [Appelbaum, 

1975].  The  self-consistent  pseudopotential  method  has  also  been 
successfully  applied  to  calculation  of  the  electronic  structure  of 
vacancy  defects  [Louie,  1976a],  metal-semiconductor  interfaces 
[Louie,  1976b],  and  even  molecules  [M . L.  Cohen,  1975]. 

C.  Pseudopotentials  in  Other  Fields 

The  use  of  the  pseudopotential  concept  is  by  no  means  limited 
to  solid-state  band  structure  calculations.  However,  to  make  con- 


tact  with  other  applications,  we  must  first  consider  a formulation 
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due  to  Heme  [1970].  It  makes  use  of  the  fact  that,  xf  we  know  the 
amplitude  for  scattering  an  electron  off  the  closed-shell  ion,  we 

have  all  of  the  information  needed  to  calculate  the  band  structure.  | 

These  amplitudes  may  be  expressed  as  the  phase  shifts  for  each 

partial  wave  as  a function  of  energy.  If,  for  a given  £,  there  are 

n.£  electron  states  occupied  in  the  ion,  the  scattering  states  must 

have  ri£  radial  nodes  in  order  to  be  orthogonal  to  the  occupied 

states.  In  the  scattering  amplitude,  this  is  just  a contribution  of 

n£jr  to  the  phase  shift.  We  can  remove  the  n^Tr  without  affecting 

the  scattering  amplitude,  and  then  construct  that  potential  which 

would  give  the  same  phase  shifts.  This  is  just  the  pseudopotential. 

With  this  in  mitid,  we  see  that  the  optical  potential  used  in 
nuclear  physics  is  really  a pseudopotential  [McCarthy,  1968, 
p,  255].  In  its  simplest  form,  the  optical  potential  is  a constant 
potential  inside  the  nuclear  radius,  with  a small  imaginary  part  to 
account  for  scattering  into  other  channels.  (Note  the  similarity  to 
the  Heine-Abarenkov  model  pseudopotential.)  This  is  a pseudo- 
potential in  the  sense  that  it  is  an  effective  potential  which 
approximates  a rather  complicated  situation. 

In  field  theory,  the  Lamb  shift  is  calculated  with  what  is 
effectively  a pseudopotential  technique  [Bjorken,  1964,  pp.  177-9]. 

The  electron- proton  scattering  amplitude  is  calculated,  including 


vacuum-polarization  corrections.  It  can  be  written 
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t,  f,  o 


m 


(constant) 


u(p) 


(1.2) 


which  can  also  be  derived  from  a potential 

3 

Ze^  4q  Zo  (constant)  6 (r) . (1.3) 

' r 3 m^ 

This  "pseudopotential"  is  used  to  calculate  the  energy  levels  of 
the  bound  state  problem  (the  hydrogen  atom) . 


D.  The  Heterojunction  Model 

The  model  by  which  we  have  studied  heterojunction  band 
lineups  is  logically  divisible  into  two  parts.  The  first  concerns  the 
problem  of  calculating  the  bulk  semiconductor  energy  band  locations 
with  respect  to  the  periodic  electrostatic  potential.  This  has  been 
done  by  a self-consistent  pseudopotential  technique.  The  second 
part  of  the  model  is  our  scheme  for  matching  the  electrostatic 
potentials  across  the  heterojunction. 

The  self-consistent  pseudopotential  technique  was  developed 
with  an  emphasis  on  preserving  information  about  the  electrostatic 
potential.  The  pseudopotential  of  the  ion  core  was  constructed 
from  the  electrostatic  potential  of  a charge  distribution  derived  from 
atomic  Hartree-Fock  calculations.  To  this  was  added  a so-called 
pseudizing  potential  which  approximated  the  effect  of  the  core  state 
orthogonality.  The  valence  electron  charge  density  was  derived 
from  calculated  pseudo  wavefunctions , and  this  was  used  to 
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calculate  an  electrostatic  potential  ana  an  exchange  interaction, 
using  a local  exchange  approximation.  The  calculations  were 
Iterated  until  the  valence  charge  density  self-consistently  repro- 
duced the  pseudopotential  from  which  it  was  derived. 

The  self-consistent  valence  charge  density  was  used  along 
with  the  ion  core  charge  distribution  to  calculate  the  electrostatic 
potential.  This  was  used  to  calculate  a reference  potential,  which 
was  chosen  to  be  the  average  of  the  electrostatic  potentials  at  the 
two  large  interstitial  vacancies  in  the  zinc  blende  structure.  The 
reason  for  this  choice  is  that  charge  density  calculations  [Walter, 
1971]  show  that  much  of  the  interstitial  space  in  a semiconductor 
crystal  is  nearly  empty.  Thus  much  of  the  interfaciai  area  of  a 
heterojunction,  along  which  the  electrostatic  potentials  must  be 
matched,  is  a region  of  low  charge  density.  We  therefore  think  that 
it  is  most  appropriate  to  define  a reference  potential  in  that  part  of 
the  crystal  with  low  charge  density.  With  these  considerations, 
our  final  prescription  for  predicting  heterojunction  band  lineups  is 
to  assume  that  the  reference  levels  line  up  at  the  heterojunction. 


CHAPTER  n 


CALCULATIONS 

A.  The  Ionic  Pseudopotential 

Since  we  are  interested  in  calculating  the  electrostatic 
potential  along  with  the  band  structure,  it  is  important  to  construct 
the  ionic  pseudopotential  in  such  a way  that  we  can  recover  the 
electrostatic  part.  This  consideration  has  not  been  important  in 
other  pseudopotential  work  [Appelbaum,  1973;  Schluter,  1975] 
where  the  ionic  pseudopotential  was  taken  to  be  of  some  convenient 
analytical  form  with  adjustable  parameters.  In  our  work  it  is 
desirable  to  start  with  a reasonable  approximation  to  the  electro- 
static potential  of  the  ion,  and  to  add  terms  to  approximate  the 
orthogonality  requirement. 

The  ionic  pseudopotential  needs  to  satisfy  three  criteria.  It 
must  be  sufficiently  weak  that  the  lowest  bound  electron  states 
correspond  to  the  valence,  rather  than  the  core,  electrons.  It  must 
be  a well-behaved  function,  with  no  singularities  in  magnitude  or 
slope.  This  allows  rapid  fall-off  of  the  Fourier  transform,  which  is 
necessary  for  the  proper  convergence  of  the  band  structure  calcula- 
tions. Finally,  it  must  be  localized  in  the  vicinity  of  the  ion.  This 
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imposes  a relationship  between  the  zero  wavevector  Fourier  compo- 
nent, which  does  not  affect  the  band  gaps  (but  which  does  affect 
the  band  energy  with  respect  to  the  electrostatic  potential)  , and  the 


nonzero  wavevector  components,  which  do.  Thus,  in  adjusting  the 
ionic  pseudopotential  to  fit  the  band  structure,  the  nonzero  wave- 
vector  components  act  to  determine  the  k = 0 component,  which  can 
therefore  be  regarded  as  a meaningful  quantity.  Consequently,  we 
must  use  an  analytic  form  for  the  ionic  pseudopotential,  with 
adjustability  provided  by  parameters  within  the  analytic  expression. 
We  cannot  adjust  individual  Fourier  components,  as  is  done  in  the 
EPM  and  in  some  OPW  methods  [see  Herman,  1967,  for  example]. 
We  begin  the  construction  of  the  ionic  pseudopotential  with 


the  electrostatic  potential  of  the  ion.  We  know  that  for  very  small 
distances  from  the  nucleus,  the  potential  must  be  of  the  form  Ze/r, 
where  Z is  the  atomic  number.  For  large  distances,  the  potential 
must  be  of  the  form  (Z-Q)e/r,  where  Q is  th  number  of  core  elec- 
trons. Now,  in  order  to  satisfy  the  smoothness  criterion  for  the 
pseudopotential,  the  Ze/r  behavior  at  the  origin  must  be  cancelled 
out  in  the  pseudizing  process.  This  is  easily  done  if  we  can 
assume  that  the  core  electrons  are  distributed  as 
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where  a is  a parameter  characterizing  the  size  of  the  ion. 
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We  must  now  consider  the  validity  of  this  form  for  the  core 
charge  distribution.  There  are  a number  of  criteria  that  could  be 
applied  to  test  whether  this  provides  a good  approximation  to  cal- 
culated ionic  charge  distributions,  but  the  one  which  was  chosen 
was  that  the  model  should  reproduce  the  Fourier  components  for  low 
wavevectors.  The  reason  for  this  choice  is  that  the  band  structure 
calculations  are  performed  in  k-space,  using  a limited  number  of 
plane  waves,  so  the  Fourier  components  of  the  potential  with  small 
k are  the  important  quantities.  It  is  also  very  easy  to  apply  this 
test,  because  the  Fourier  components  of  atomic  and  ionic  electron 
distributions  ("form  factors")  from  Hartree-Fock  calculations  are 
tabulated  in  the  International  Tabl  ;s  for  X-ray  Crystallography 
[MacGillayry , 1962,  pp.  201-6],  The  Fourier  transform  of  our 
assumed  charge  distribution  is 


^core^*^^ 


« Q(  1- 


\ 


Q 


(II . 2) 


which  would  fit  the  form  factors  for  just  about  any  localized  distri- 
bution, for  sufficiently  small  k.  However,  this  expression  is 
capable  of  fitting  the  tabulated  form  factors  to  within  5 percent  for 
wavevectors  less  than  four  or  five  reciprocal  Bohr  radii.  The 
largest  wavevectors  included  in  the  band  structure  calculations 
were  about  2.5  Sq  Therefore,  it  appears  that  for  our  purposes, 
the  expression  (II.  1)  for  the  ionic  charge  distribution  is  quite  ade- 
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The  pdrameier  u in  (II . 1)  must  be  determined  from  Hartree- 
Fock;  calculations  for  the  ions.  The  International  Tables  for  X-rny 
Crystallography  do  not  list  form  factors  for  all  of  the  elements  of 
interest  in  semiconductor  work,  so,  having  demonstrated  the  validity 
of  our  assumed  form,  it  was  necessary  to  find  a more  complete  set 
of  data.  This  is  readily  provided  by  the  work  of  Desclaux,  where 
the  results  of  relativistic  Hartree-Fock  calculations  are  presented 
for  atoms  with  Z=1  to  Z=120  [Desclaux,  1973].  It  is  most  impor- 
tant to  do  all  of  our  calculations  in  exactly  the  same  way,  so  we 
have  taken  aJT  of  our  ionic  data  from  this  work.  The  charge  distri- 
bution data  in  Desclaux  are  presented  in  terms  of  various  radial 
moments  for  each  atomic  shell.  We  used  the  average  second  moment 
(r^)  for  the  occupied  shells  of  the  ion  to  determine  a according  to 
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(II.  3) 


The  values  of  a thus  determined,  along  with  the  other  parameters  of 
the  ionic  pseudopotential,  are  given  in  Table  I. 

The  electrostatic  potential  energy  of  an  electron  in  the  field 
of  the  ion,  using  our  model  charge  distribution,  is  (in  c.g.s.  nota- 
tion) 


Vpjr)  = -7- 
es  r 


(Z-Q)  + Qe  . 


(II.  4) 


As  mentioned  earlier,  the  Ze/r  singularity  at  the  origin  must  be 
removed  in  the  pseudizing  process.  This  is  easily  done  by  choosing 
part  of  the  pseudizing  potential  to  be 
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TABLE  I 

IONIC  PSEUDOPOTENTIAL  PARAMETERS 


Element 

Z 

Q 

a (a.u.) 

7 (a.u.) 

Vq  (Ry.) 

A1 

13 

10 

4.04 

1.7 

43 

Si 

14 

10 

4.52 

2.25 

44 

P 

15 

10 

4.99 

2.5 

41 

S 

16 

10 

5.47 

2.6 

41 

Zn 

30 

28 

3.33 

2.0 

25 

Ga 

31 

28 

3.64 

2.0 

28 

Ge 

32 

28 

3.94 

2.06 

36.8 

As 

33 

28 

4.22 

2.2 

42 

Se 

34 

28 

4.49 

2.25 

44 

Cd 

48 

46 

3.02 

1.8 

46 

In 

49 

46 

3.24 

2.0 

52 

Sb 

51 

46 

3.63 

1.8 

56 

Te 

52 

46 

3.81 

1.8 

58 
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The  sum  of  V (r)  and  V , (r)  is  a potential  well  that  is  smoothly 

0 5 JD  S 1 

2 

rounded  at  the  origir»,  as  shown  in  Figure  2,  with  a -(Z-O)e  /r 
behavior  at  large  r,  as  desired.  There  are  no  adjustable  parameters 
in  Vpgj , and  it  leaves  the  potential  at  the  origin  still  fairly  strongly 
attractive,  so  we  must  add  another  term  to  produce  a satisfactory 
ionic  pseudopotential.  This  term  was  chosen  to  be  of  the  form 


^ps2^''^  ^ [2^^ 
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(II.  6) 


where  Vg  and  7 are  freely  adjustable  parameters.  The  complete 
ionic  pseudopotential  for  Ge  is  also  shown  in  Figure  2.  The  Fourier 
transform  of  the  ionic  pseudopotential  is 


V.  (k)  = 
ion 


es 


V , (k) 
psl 


_ 87t  r (Z-Q)  Q 1 ^ 87t  Z ^ -(kV2y^) 

k^+Q^J  « kV(Q/Z)a^  n 

(II.  7) 

where  is  the  unit  cell  volume  (the  appropriate  normalization  in  a 
periodic  crystal)  and  the  Stt  really  represents  Aire^ , but  e^=2  when 
the  chosen  units  are  Rydbergs  and  Bohr  radii. 

The  expression  (II. 7)  diverges  as  k-0,  as  any  potential 
arising  from  a nonzero  net  charge  must.  In  the  semiconductor  crys- 
tal, this  poses  no  real  problem,  since  the  Coulomb  potential  is 


Figure  2.  The  components  of  the  ionic  pseudopotential  for  Ge 
The  lowest  curve  is  the  electrostatic  potential;  the 
center  curve  is  the  electrostatic  potential  after  the 
singularity  has  been  removed  by  the  addition  of  Vp^ 
The  upper  curve  is  the  final  pseudopotential. 
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screened  by  the  valence  electrons,  but  we  need  a well-defined, 
finite  value  for  for  the  band  structure  calculations.  The 

solution  is  to  introduce  a fictitious  screening  charge  distribution  in 
the  V^Qj^(k=0)  calculation  for  bookkeeping  purposes,  and  to  sub- 
tract a corresponding  charge  distribution  from  the  valence  band 
charge.  This  charge  distribution  was  chosen  to  be  of  the  form 


Pp(r) 


(Z-0) 


(2;r)V^ 


(II.  8) 


where  (i  is  an  arbitrary  parameter.  When  the  effect  of  this  ficti- 
tious charge  distribution  is  taken  into  account,  the  (Z-Q)A^  term 
in  (II.  7)  IS  replaced,  for  k=0 , by  (Z-0)/2^3^.  This  scheme  will  be 
discussed  in  greater  detail  when  we  consider  the  actual  calculation 
of  the  electrostatic  potential. 

The  initial  values  of  the  parameters  Vq  and  y were  chosen  by 
fitting  the  energy  levels  of  the  closed-shell-plus-one-electron  ion. 
The  parameters  were  then  adjusted  in  response  to  the  band  structure 
results,  until  an  acceptable  fit  was  obtained.  The  final  ionic 
pseudopotential  was  again  used  to  calculate  the  ionic  spectrum,  to 
ensure  that  the  pseudopotential  was  a reasonable  representation  of 
the  ion.  The  method  of  calculating  the  ionic  spectra,  and  the  final 
results,  are  given  in  Appendix  A. 
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B.  Self-Consistem  Band  Structure  Cdlculdtions 

With  our  ionic  pseudopotential,  we  can  calculate  the  energy 
band  structure  self-consistently . The  procedure  is  to  choose  an 
initial  valence  charge  distribution,  derive  from  it  the  electrostatic 
and  exchange  potentials,  and  add  these  to  the  ionic  pseudopotential 
to  obtain  the  total  pseudopotential.  This  is  then  used  to  calculate 
the  wavefunctions  for  a suitable  number  of  states.  The  appropriate 
linear  combinations  of  the  absolute  squares  of  the  wavefunctions 
give  the  valence  electron  distribution.  This  is  used  to  generate  a 
new  pseudopotential,  and  the  entire  process  is  repeated  until  the 
pseudopotential  converges  to  a consistent  value. 

The  choice  of  the  initial  valence  electron  distribution  is  not 
particularly  critical.  Appelbaum  and  Hamann  [1973]  simply  assumed 
a linear  dielectric  response  to  the  ionic  pseudopotential.  We  have 
tried  to  speed  convergence  by  using  an  initial  valence  charge 
density  calculated  from  published  empirical  pseudopotentials.  In 
certain  cases,  where  such  potentials  were  not  readily  available, 
we  have  started  the  calculation  with  the  self-consistent  charge 
distribution  for  a related  compound.  (For  example,  we  initiated  the 
CdSe  calculation  with  the  ZnSe  charge  distribution.) 

Since  the  charge  distribution  is  represented  as  a Fourier 


series,  calculating  the  electrostatic  (Hartree)  potential  is  a 
straightforward  matter  of  dividing  each  component  by  the  square  of 
its  wavevector.  The  exchange  interaction  is  more  complicated. 
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We  have  usee  the  generalized  Slater  approximation,  commonly 
known  as  the  "Xu"  approximation  [Slater,  1974].  In  this  approach, 
the  exchange  interaction  is  approximated  by  a local  potential  which 


IS  proportional  to  the  cube  root  of  the  charge  density.  In  our  cal- 
culations, the  exact  evaluation  of  this  approximation  would  require 
summing  the  Fourier  series  for  the  charge  density  at  a number  of 
points,  taking  the  cube  root  at  each  point,  and  re-analyzing  the 
series.  Instead  of  this  rather  time-consuming  procedure,  we  have 
used  a polynomial  approximation  to  the  cube  root,  which  can  be 
evaluated  by  direct  manipulations  on  the  series.  The  valence 
charge  densities  in  semiconductors  are  limited  to  a range  of  about 
one  to  forty  electrons  per  unit  cell  [Walter,  1971].  Fitting  a quad- 
ratic polynomial  to  the  cube  root  over  this  range  gives 

= 0.847  + 0.142P  - 0.00209P  (II. 9) 

which  IS  shown  in  Figure  3.  In  the  first  round  of  calculations  we 
used  a pure  Slater  approximation  with  a=l,  but  for  the  final  calcu- 
lations we  used  a=0.8.  This  value  has  been  used  successfully  in 
surface  calculations  [Schlliter,  1975],  and  it  gave  improved  band 
structures  in  our  own  calculations. 

The  pseudo  wavefunctions  and  their  energies  were  calculated 


using  established  pseudopotential  techniques.  The  Hamiltonian 
matrix  was  set  up  on  a basis  of  plane  waves  Ik+S)  where  k is  a 
vector  in  the  reduced  Brillouin  zone  and  S is  a reciprocal  lattice 
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Figure  3.  The  cube  root  of  the  charge  density  (solid  line)  and  the 
quadratic  approximation  (11.9)  (dotted  line),  and  the 


range  over  which  this  approximation  was  actually  used. 
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vector.  For  each  value  of  k,  all  plane  waves  such  that  lk+GI^=E^ 
were  included  m the  direct  calculation,  and  the  remaining  plane 
waves  such  that  lk+5l^  = E2  were  included  via  the  second-order 
perturbation  scheme  of  Lowdin  [1951]  as  modified  by  Brust  [1964]. 
We  found  that  setting  E^=7  and  ££  = 14  (in  units  of  (27r/a)^)  yielded 
good  band  structure  results  without  requiring  excessive  computa- 
tion. The  criterion  for  "good"  band  structures,  in  this  regard,  is 
that  the  band  energies  should  not  change  excessively  as  Ic  moves 
through  a point  where  the  set  of  included  G's  changes.  If  the 
basis  set  is  too  small,  the  sudden  inclusion  of  different  pseudo- 
potential  components  leads  to  significant  changes  in  the  eigen- 
values. The  eigenvalues  were  calculated  by  numerically 
diagonalizing  the  Hamiltonian,  and  the  wavefunctions  were  derived 
from  the  eigenvectors  of  the  foldec -down  Hamiltonian  by  a partial 
unfolding  of  the  Lowdin  scheme. 

The  pseudocharge  density  was  calculated  using  the  Ghadi- 
Cohen  three-point  Brillouin  zone  sample  [Chadi,  1973].  This  has 
been  shown  to  give  charge  densities  which  agree  to  within  about 
1 percent  with  calculations  involving  many  more  points  in  the 
Brillouin  zone.  The  expression  used  to  calculate  the  charge  density 
is 

P(r)  = + ^p^(F) 


(II.  10) 
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where  the  are  the  charge  densities  from  the  various  points, 

symmetrized  with  respect  to  the  point-group  operations: 

V'  1 V'  ( i)  ■ 2 

= L TiT  L (II.  11) 

where  the  index  i refers  to  the  band  and  is  summed  over  the 
valence  bands,  R is  an  element  of  the  crystal  point  group  G which 
has  nQ  such  elements.  The  k vectors  used  are 

A = ^ (I-,  0 , 0) 

a 

W = ^ (1,  I , 0)  (11.12) 

a 

V — - n\ 

“ - g i2  / 2 I u; . 

The  simple  procedure  of  iterating  the  band  structure  and 
charge  density  calculations  converges  quickly  for  the  Group  IV 
semiconductors,  but  for  the  heteropolar  (III-V  and  II-VI)  semi- 
conductors the  convergence  is  oscillatory.  This  was  easily  damped 
out  by  taking  a linear  combination  of  the  "new"  and  "old"  charge 
densities  as  the  input  to  the  next  iteration;  a contribution  of  one- 
half  from  each  worked  very  well.  The  average  root-mean-square 
change  in  the  electrostatic  potential  was  calculated  at  each  itera- 
tion. The  iterations  were  terminated  when  this  reached  a level  of  | 

] 
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0.001  or  0.002  Rydbergs. 
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C.  Evaluation  of  the  Electrostatic  Potential 

With  the  ion  model  and  the  valence  charge  distribution  from 
the  band  structure  calculations,  we  have  complete  information  on 
the  microscopic  charge  distribution  in  the  crystal.  In  principle  the 
calculation  of  the  electrostatic  potential  is  a straightforward  solu- 
tion of  Poisson's  equation.  In  practice,  some  care  must  be  taken 
if  the  calculation  is  to  produce  meaningful  results.  The  difficulties 
arise  in  the  methods  of  representing  the  potential,  so  we  must 
digress  briefly  to  discuss  these  problems. 

There  are  fundamentally  two  ways  to  represent  a periodic 
scalar  function  (such  as  the  charge  density  or  the  potential)  in  a 
crystal  [Kittel,  1963,  ch.  1].  A function  F can  be  represented  as  a 
Fourier  series  in  the  reciprocal  lattice  vectors  G: 

F(T)  = a.e^^i'"'  . (II.  13) 

i 

The  same  function  can  also  be  represented  as  a sum  of  "atomic" 
functions  f centered  at  lattice  points  R^: 

F(r)  = Z (II.  14) 

i 


The  relationship  between  the  two  representations  is 
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^i  ~ n 1*  ' • (II-  15) 

In  applying  these  representations  to  numerical  computations. 
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we  see  that  each  is  suited  to  a particular  kind  of  function.  The 
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Fourier  series  works  well  for  smooth  functions  which  need  not  be 
well  localized.  However,  it  is  incapable  of  representing  singulari- 
ties when  truncated  at  any  finite  number  of  terms,  as  it  must  be  in 
any  practical  computation.  The  atomic  functions,  on  the  other 
hand,  may  be  quite  singular,  but  they  must  be  well  localized. 

In  the  case  of  the  crystal  electrostatic  potential,  we  make 
use  of  both  representations.  The  Coulomb  singularity  at  the 

< 

nucleus  cannot  be  treated  with  a truncated  Fourier  series.  On  the 
other  hand,  the  valence  charge  distribution  exists  as  a Fourier 
series  and  representing  it  in  terms  of  atomic  functions  would 
require  some  very  tedious  fitting  to  Equation  (11.15).  Therefore,  we 
use  atomic  functions  to  represent  the  ionic  charge  distribution, 
while  a Fourier  series  is  used  for  the  valence  charge.  This  creates 

a problem  in  itself,  for  we  cannot  solve  Poisson's  equation  in  ■ 

either  representation  if  the  net  charge  involved  is  nonzero.  With 
the  Fourier  series,  this  means  that  the  G=0  component  is  nonzero. 

Dividing  by  G gives  an  undefined  result.  In  the  atomic-function 
representation,  a net  charge  implies  a Coulomb  tail.  Since  the 
number  of  lattice  points  at  a radius  R is  proportional  to  R , the 
1/R  Coulomb  term  is  not  sufficient  to  cause  the  series  to  converge. 

These  difficulties  are  clearly  an  artifact  of  the  mathematical  formu- 
lation, because  the  crystal  itself  is  neutral.  In  fact,  they  are  a 
corollary  to  the  problem  of  the  divergence  of  V^qj^(K=0)  discussed 


earlier. 
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The  soluuon  to  these  difiiculties  is  to  remove  en  arbitrary 
charge  distribution  from  the  valence  charge  and  add  it  onto  the  icn, 
so  that  both  the  valence  and  ionic  distributions  are  neutral  overall. 
This  fictitious  distribution  is  given  in  Equation  (II.  8).  (The  choice 
of  a Gaussian  is  essentially  required  because  we  need  rapid  con- 
vergence in  both  r and  k space.)  The  Fourier  transform  of  this 
charge  distribution  (with  the  appropriate  structure  factors)  is  sub- 
tracted from  the  valence  charge,  allowing  us  to  solve  Poisson's 
equation  and  set  the  G=0  Fourier  component  of  the  potential  equal 
to  zero  (the  usual  boundary  condition)  . The  fictitious  charge  dis- 
tribution is  added  to  the  ionic  charge  distribution,  cancelling  the 
Coulomb  tail  and  allowing  us  to  calculate  the  ionic  contribution  to 
the  electrostatic  potential  according  to  (11.14).  The  two  contribu- 
tions are  then  added  to  obtain  the  total  electrostatic  potential. 

There  is  one  remaining  subtle  point  in  this  procedure.  The 
boundary  conditions  are  defined  so  that  the  Fourier  series  makes  no 
contribution  to  the  average  (C=Q)  electrostatic  potential.  However, 
the  fictitious  charge  plus  the  (Z-Q)  part  of  the  nuclear  charge  con- 
tribute a term 


^ IZzgj 

Q 2/3^ 


(II.  16) 


to  the  average  electrostatic  potential,  which  must  be  included  in 
the  G=0  component  of  the  pseudopotential,  as  noted  previously. 
Thus  the  arbitrary  parameter  fi  affects  both  the  electrostatic 
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potential  and  the  band  energy,  but  it  cannot  affect  the  relation 
between  the  two.  (It  effectively  sets  the  zero  on  the  energy  scale.) 
The  existence  of  this  arbitrary  but  irrelevant  parameter  has  proved 
to  be  quite  useful  in  detecting  program  errors.  (If  it  is  not  irrele- 
vant, there  is  an  error.) 

There  are  practical  limits  on  the  value  of  imposed  by  the 
desire  for  rapid  convergence.  In  the  Fourier  series,  we  included 
G's  up  to  about  3 a.u.  If  we  want  the  terms  at  the  cutoff  to  be 
smaller  than  10"“^  Ry.  , we  should  use  /3<1.8.  In  the  summation 
over  a cluster  of  lattice  points,  the  cutoff  was  at  a distance  of 
10  a.u.  Applying  the  same  criterion  for  the  size  of  the  terms  at 
cutoff,  we  find  that  we  should  use  /3>0.3.  Most  of  the  calcula- 
tions actually  used  ji=0.5  a.u. 

D.  Matching  Scheme 

The  matching  scheme,  which  specifies  how  the  electrostatic 
potentials  should  line  up  at  the  heterojunction,  could  be  developed 
at  a number  of  levels  of  sophistication.  In  the  earlier  stages  of 
this  investigation,  we  thought  that  it  would  probably  be  necessary 
to  model  the  detailed  charge  distribution  at  the  interface,  probably 
by  assuming  some  sort  of  dielectric  response.  The  extreme  case  of 
a matching  scheme  would  be  to  do  the  full  quantum-mechanical 
interface  calculation,  as  is  currently  being  done  for  surfaces 
[Appelbaum,  1975;  Schlliter,  1975].  We  have  found , however,  that 


29 

a very  simple  matching  scheme  gives  remarkably  good  heterojunc- 
tion  predictions. 

Our  matching  scheme  is  suggested  by  considering  a muffin- 
tin  potential  [Harrison,  1970,  p.  87].  This  potential,  which  is 
commonly  used  for  APW  band  structure  calculations,  consists  of 
spherically  symmetric  wells  within  a certain  core  radius  from  the 
ions,  and  is  constant  everywhere  else.  If  the  electrostatic  poten- 
tial in  a semiconductor  were  really  of  a muffin-tin  form,  then  the 
matching  scheme  would  be  completely  obvious;  the  flat  parts  of 
the  potentials  would  line  up. 

Our  calculations  indicate  that,  for  the  cubic  semiconductors, 
the  electrostatic  potential  is  indeed  flat  to  within  about  a volt 
inside  the  large  interstitial  spaces  in  the  crystal  structure.  These 
spaces  are  centered  on  the  points  a(LLi)  and  a(|,|,|)  if  the 
atoms  are  at  a (0 , 0 , 0)  anda(i,i,j).  For  the  zincblende  semi- 
conductors the  two  interstitial  points  are  not  equivalent,  but  the 
potential  difference  is  typically  less  than  one  volt.  Therefore,  we 
have  simply  taken  the  average  of  the  potentials  at  the  two  inter- 
stitial points  as  our  reference  potential.  Then,  following  the 
muffin-tin  analogy,  the  matching  scheme  is  simply  to  line  up  the 
reference  potentials. 

This  scheme  can  be  expected  to  work  best  for  heterojunctions 
grown  on  nonpolar  {110}  crystal  faces.  In  this  case,  the  interstitial 
points  lie  in  the  atomic  planes  adjacent  to  the  junction,  with  equal 
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numbers  of  each  type.  For  a heterojunction  on  a {111)  polar  face, 
the  interstitial  points  are  slightly  displaced  from  the  heterojunction 
plane,  with  the  two  types  displaced  in  opposite  directions.  There- 
fore, it  would  probably  be  more  appropriate  to  use  a slightly 
different  weighting  factor  for  the  two  points  in  such  cases.  Experi- 
mentally, the  only  investigation  of  orientation  effects  is  that  by 
Fang  and  Howard  [1964].  They  found  differences  of  about  0.2  eV  in 
band  lineup  between  the  two  (ill)  orientations  of  Ge-GaAs  junc- 
tions. This  probably  indicates  the  accuracy  we  can  expect  from  our 
matching  scheme  as  long  as  we  neglect  orientation  effects. 

It  should  be  emphasized  that  this  matching  scheme  is  pri- 
marily appropriate  to  ideally  abrupt  heterojunctions.  If  the  compo- 
sition is  graded  over  more  than  a few  atomic  layers  , the  reference 
levels  of  the  pure  semiconductors  cease  to  be  important  and  the 
equalization  of  the  Fermi  levels  begins  to  dominate.  Also,  this 
matching  scheme  explicitly  neglects  interface  states.  We  believe, 
for  physical  reasons,  that  the  interface  states  at  heterojunctions 
can  result  only  from  imperfections  and  dislocations  at  the  junction. 
Therefore,  our  model  is  limited  to  predictions  for  well  lattice- 
matched  heterojunctions , 


CHAPTER  III 
RESULTS 

A.  Band  Structure  Results 

We  have  applied  the  model  described  in  the  previous  chapter 
to  fifteen  semiconductors  with  cubic  structures.  The  ionic  pseudo- 
potentials were  adjusted  to  fit  the  band  structures,  with  particular 
emphasis  on  the  fundamental  direct  and  indirect  gaps,  but  the 
pseudopotential  for  a given  ion  was  the  same  for  all  compounds 
involving  that  ion.  Of  the  Group  IV  semiconductors,  only  silicon 
and  germanium  were  calculated,  since  these  are  the  ones  of  tech- 
nological interest.  There  is  a great  deal  of  interest  in  the  III-V 
quaternary  alloy  systems  (Ga  , In)  (As  , P)  and  (A1 , Ga)  (As  , Sb) , so  we 
calculated  the  pure  compounds  for  each  of  these  systems.  Of  the 
II-VI  compounds,  we  calculated  all  of  the  possible  combinations  of 
Zn  and  Cd  with  S,  Se,  and  Te. 

Table  II  gives  the  band  gaps  which  are  usually  considered  in 
discussing  band  structures.  The  experimental  data  are  given  for 
comparison  with  the  calculated  values.  The  symmetry  labels  are 
those  appropriate  to  the  zincblende  structure.  In  quoting  the 
experimental  data,  we  have  chosen  electroreflectance  data  when 
available,  preferably  at  low  temperature.  We  have  used  the 
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"traditional"  assignment  of  the  optical  structures.  The  Eq  and  Eq 
peaks  in  Cardona's  notation  [Cardona,  1969J  are  assigned  to  the 
Tisv-r^c  and  r'i5v“’^15c  transitions,  respectively.  The  E^  and 
Ej^'  sets  of  peaks  are  assigned  to  and  A^^'A^^  transi- 
tions, and  the  E2  peaks  are  assigned  to  transitions  from  to 

Xi-  and  Xo  . With  the  refinement  of  both  theory  and  experiment, 

1C  oc 

it  now  appears  that  this  scheme  is  too  naive.  The  optical  structures 
arise  from  transitions  over  sizeable  volumes  of  the  Brillouin  zone, 
and  detailed  critical-point  analyses  are  needed  to  properly  interpret 
the  optical  structure  [see  Saravia  , 1968,  for  example].  However, 
the  assigned  locations  in  the  Brillouin  zone  have  not  changed  radi- 
cally, and  we  think  that  the  older  interpretation  is  adequate  for  our 
work . 

We  made  no  provision  for  calculating  the  spin-orbit  splitting , 
so  the  experimental  data  have  been  corrected  for  this.  To  all  gaps 
involving  ^^5^/  one-third  of  the  spin-orbit  splitting  of  this  level 
(Aq)  has  been  added.  The  A 3^  levels  have  been  corrected  by 
adding  one-half  of  A^.  The  spin-orbit  splitting  at  X is  almost 
always  negligible. 

The  calculated  values  given  in  Table  II  have  been  obtained 
from  calculations  sampling  the  Brillouin  zone  primarily  along  the 
A (111)  and  A ( 100)  directions . The  r point  is  the  center  of  the 
zone.  The  indirect  gaps,  and  are  from  the  top 

of  the  valence  band  to  the  local  minima  along  the  respective 
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directions.  The  direct  gaps  along  A were  obtained  by  finding  either 
the  minimum  gap  or  a region  where  the  gap  is  approximately  con- 
stant, in  an  attempt  to  guess  where  the  critical  point  might  be. 
Accordingly,  these  gaps  are  only  given  to  0. 1 eV.  In  the  A direc- 
tion, the  direct  gaps  are  calculated  at  the  point  X. 

An  examination  of  Table  II  shows  that  we  have  generally  been 
able  to  fit  the  first  direct  gap  at  r and  the  lowest  gap  of  the  indi- 
rect gap  materials  to  within  about  0.1  eV.  The  worst  fit  to  these 
levels  is  that  for  AlAs,  where  the  direct  gap  is  0.2  eV  low  and  the 
indirect  gap  is  0.3  eV  low.  (However,  the  experimental  data  were 
taken  on  n-type  samples  with  carrier  concentration  in  the  high- 
10^®-cm~^  range  [Monemar,  1973],  so  there  could  possibly  be  a 
Burstein  shift  of  the  order  of  0.1  eV  in  this  gap.)  Other  indirect 
gaps  show  a reasonable  agreement  with  the  data,  except  the 
Fisv-Aic  gap  in  InAs.  The  data  for  InAs  are  from  high-pressure 
Hall-effect  experiments,  which  found  a local  minimum  in  the  con- 
duction band  0.7  eV  above  the  band  edge,  which  was  interpreted  to 
be  along  the  A direction.  However,  it  is  difficult  to  see  how  such 
a minimum  might  occur,  considering  the  direct  gaps  and  typical 
curvatures  of  the  upper  valence  band,  so  we  should  not  be  too 
worried  by  this  discrepancy. 

The  gaps  generally  fit  the  experimental  data  to  about 

0.3  eV  and  the  same  may  be  said  of  the  gap.  although 


there  are  a few  cases  where  the  disagreement  is  considerably 
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larger.  The  upper  conduction  bands  are  not  very  well  fitted;  the 

average  error  for  the  A_  -A.  gap  is  about  0.5  eV.  The  direct  gaps 

iv  3c 

at  X are  particularly  difficult  to  fit  with  our  model.  The  disagree- 
ment in  the  -X,  gap  runs  around  1 eV  and  there  are  some  par- 
5v  Ic 

ticularly  glaring  errors  in  the  gap  between  the  two  conduction  bands 
at  X. 

It  appears  that  no  experimental  data  exist  for  cubic  cadmium 
selenide,  but  the  band  gap  is  in  reasonable  agreement  with  that  of 
the  hexagonal  crystal.  For  those  II-VI  compounds  which  exist  in 
both  structures,  the  band  gaps  are  usually  about  the  same. 

The  remarkable  thing  about  our  model  is  that  we  can  fit  the 
fundamental  gaps  for  this  number  of  semiconductors  with  the  same 
ionic  pseudopotentials  (which  have  only  two  adjustable  parameters). 

There  are  limits  to  this,  however.  We  have  the  pseudopotentials 
necessary  to  calculate  AlP  and  InSb,  but  these  give  direct  band  gaps 
which  are  in  error  by  0.6  and  0.8  eV,  respectively. 

The  object  of  our  band  structure  calculations  has  of  course 
been  to  calculate  the  energies  of  the  bands  in  relation  to  our  refer- 
ence potential.  These  results  are  given  in  Table  III.  In  this  case  I 

I 

the  calculations  have  been  corrected  for  spin-orbit  splitting  by  | 

adding  one-third  of  the  experimentally  observed  splitting  Aq  to  the 

. 1 

valence  band  energies.  Also  given  is  the  difference  in  the  electro-  ! 

1 

static  potential  (AV)  between  the  two  interstitial  points  used  to  i 

define 

I 

i 

^ 
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TABLE  III 

BAND  EDGE  ENERGIES  WITH  RESPECT  TO  THE 
REFERENCE  POTENTIAL 


The  energies  of  the  edge  of  the  valence  band,  including  spin-orbit 
splitting,  and  the  conduction  bands  are  given.  AV  is  the  difference 
in  potential  between  the  two  interstitial  points. 


Semiconductor 

E 

V 

E^  (direct) 

E^(indirect) 

AV 

Si 

-3.16 

-0.21 

-2.10 

0 

Ge 

-3.25 

-2.39 

-2.51 

0 

AlAs 

-3.96 

-1.21 

-2.00 

0.74 

AlSb 

-3.94 

-1.66 

-2.35 

0.42 

GaP 

-4.12 

-1.26 

-1.92 

0.66 

GaAs 

-3.96 

-2.44 

— 

0.48 

GaSb 

-3.98 

-3.04 

— 

0.19 

InP 

-4.58 

-3. 15 

— 

0.81 

InAs 

-4.38 

CO 

CO 

CO 

-- 

0. 64 

ZnS 

-5.34 

-1.53 

— 

0.92 

ZnSe 

-5.07 

-2.33 

-- 

0.74 

ZnTe 

-4.74 

-2.51 

-- 

0.46 

CdS 

-5.42 

-3.  11 

— 

1.07 

CdSe 

-5.29 

-3.53 

_ — 

0.99 

CdTe 


-4.90 


-3.41 


0.71 
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B.  Heterojunction  Predictions 

For  purposes  of  making  heterojunctions , the  semiconductors 
we  have  considered  can  be  divided  into  four  groups,  based  on  lat- 
tice constant.  Those  semiconductors  with  a lattice  constant  of 


o 

about  5.45  A are  Si,  GaP,  and  ZnS. 


o 

Around  5.65  A we  have  Ge, 


O 


GaAs , AlAs , and  ZnSe.  There  are  only  InP  and  CdS  at  5.85  A,  but 
the  group  around  6.  1 A includes  GaSb,  AlSb,  InAs , ZnTe , and  CdSe. 
If  a heterojunction  between  pure  semiconductors  is  to  be  lattice- 
matched,  both  semiconductors  must  be  chosen  from  the  same 


group. 

The  relative  band  edge  energies  for  the  silicon  group  are  shown 
in  Figure  4.  Our  model  predicts  that  most  of  the  discontinuity 
will  be  in  the  valence  band  for  all  heterojunctions  within  this  group. 
The  best  experimental  data  available  for  this  group  are  those  of 
Zeidenbergs  and  Anderson  [1967],  who  fabricated  isotype  (n-n) 
heterojunctions  between  Si  and  GaP.  These  junctions  exhibited 
ohmic  current-voltage  characteristics  at  temperatures  as  low  as 
77  K.  This  indicates  that  there  is  a negligible  conduction  band 
discontinuity,  in  qualitative  agreement  with  our  predictions. 

Lendvay,  et  al.,  [1970]  have  prepared  n-n  Si-ZnS  heterojunctions, 
and  measured  their  electrical  properties.  These  junctions  do  show 
a diode  characteristic,  with  a polarity  such  that  the  ZnS  conduction 
band  must  lie  above  that  of  Si,  but  the  data  are  not  sufficient  to 
determine  the  height  of  the  barrier. 
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The  greatest  amount  of  experimental  worx  has  been  done  on 
heterojunctions  within  the  Ge  group.  The  Ge-GaAs  system  has 
been  one  of  the  most  widely  studied  heterojunctions.  Unfortunately, 
the  data  scatter  widely.  The  more  recent  data,  however,  tend  to 
indicate  a small  conduction  band  discontinuity.  Riben  and  Feucht 
get  a value  of  0.11  eV  [Riben,  1966],  in  good  agreement  with  our 
value  of  0.07  eV  (Figure  5). 

ZnSe  has  been  used  as  an  emitter  on  both  Ge  and  GaAs  tran- 
sistors [Sieger,  1970]  but  these  devices  do  not  really  tell  us  much 
about  the  band  lineup.  Mach,  et  al,,  [1970]  have  studied  n-n 
GaAs-ZnSe  junctions  which  show  nonohmic  behavior.  In  these 
devices  a greater  current  is  passed  if  the  ZnSe  is  made  positive, 
indicating  that  the  ZnSe  conduction  band  is  lower  than  that  of  GaAs. 
However,  the  current-voltage  characteristics  are  far  from  the  expo- 
nential expected  for  a simple  diode,  so  these  data  cannot  be  con- 
sidered conclusive. 

Undoubtedly  the  best  experimental  data  available  are  those  of 
Dingle,  et  al.,  [1974  and  1975]  for  the  GaAs-AlAs  system.  They 
fabricated  periodic  arrays  of  heterojunctions  between  layers  suffi- 
ciently thin  that  they  exhibited  quantized  "square  well"  electron 
and  hole  states.  These  states  were  observed  spectroscopically  and 
from  their  energies,  the  barrier  heights  (band  discontinuities)  were 
determined.  These  structures  were  grown  by  molecular  beam 
epitaxy  and  are  very  abrupt;  the  grading  distance  is  less  than  5 A. 
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These  heterostructares  were  actually  fabricated  of  GaAs  anc 
Alj^Ga^-j^s  with  x in  the  range  0.2  to  0.3.  The  observed  valence 
band  discontinuity  is  0.15±0.03  of  the  total  difference  in  band- 
gaps.  If  we  assume  that  tne  band  energies  interpolate  linearly,  this 
ratio  should  be  constant,  independent  of  composition.  Our  calcu- 
lations give  zero  valence  band  discontinuity  between  pure  GaAs  and 
AlAs . While  this  prediction  is  not  quite  so  spectacularly  in  agree- 
ment with  the  data  as  that  which  we  obtained  in  earlier  calculations 
[Frensley,  197  6],  the  agreement  is  still  quite  reasonable. 

For  the  InP-GdS  heterojunction  we  predict  almost  continuous 
conduction  bands  (Figure  6).  The  only  experimental  data  available 
are  those  of  Shay,  et  al.,  [1976].  They  have  performed  capacitance- 
voltage  measurements  on  p-n  junctions,  and  the  results  indicate 
that  the  InP  conduction  band  lies  0.56  eV  higher  than  the  GdS  con- 
duction band.  They  further  support  this  result  by  making  theoreti- 
cal estimates  of  the  electron  affinities.  The  heterojunctions  were 
grown  by  vacuum  evaporation  at  quite  low  substrate  temperatures 
(200-250*C),  so  capacitance-voltage  measurements  should  be 
reliable.  The  origin  of  the  discrepancy  is  possibly  the  presence  of 
an  oxidized  layer  at  the  interface.  The  substrates  were  exposed  to 
the  atmosphere  while  being  prepared,  so  such  a layer  would  have 
formed.  Kroemer  [1975]  has  noted  that  the  one  case  in  which  the 
electron  affinity  rule  can  be  expected  to  hold  is  when  an  oxidized 
layer  is  present.  Another  possibility  is  that  our  calculations  may 


t 
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not  be  accurate  m this  case.  The  CdS  band  structure  is  one  of  our 
poorer  fits  to  the  experimental  data,  and  the  pseudopotential  calcu- 
lation for  the  Cd  ionic  spectrum  (Appendix  A)  , although  an  otnerwise 
excellent  fit,  shows  a large  error  for  the  5s  state.  The  InP-CdS 
heterojunction  is  something  of  an  extreme  case,  since  it  involves  a 
very  ionic  wurtzite  semiconductor.  Therefore,  it  should  not  be  too 
surprising  that  our  model  is  apparently  not  very  successful  in  this 
case . 

Some  of  the  most  interesting  band  lineups  occur  between  semi- 
conductors  of  the  6.  1 A lattice  constant  group  (Figure  7).  However, 
there  is  very  little  information  available  concerning  these  lineups. 
Much  of  the  experimental  data  have  been  interpreted  in  terms  of 
band  lineups  derived  from  the  electron  affinity  rule.  One  exception 
is  the  work  of  Nakai,  et  al.,  [1976].  They  fabricated  p-p  ZnTe-GaSb 
heterojunctions  and  studied  photoresponse  as  a function  of  photon 
energy.  They  interpreted  their  results  in  terms  of  a graded-gap 
heterojunction  with  an  overall  valence  band  barrier  of  1.0  eV.  This 
compares  favorably  with  our  predicted  value  of  0.7  6 eV,  but  one 
should  remember  that  our  model  is  not  particularly  suited  to  graded 
junctions. 

C . Heterojunctions  Involving  Semiconductor  Alloys 


There  is  a great  deal  of  current  interest  in  heterojunctions 
involving  semiconductor  alloys.  The  reason  for  this  is  that  the 
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composition  of  the  alloy  can  be  adjusted  so  as  to  achieve  a very 
good  lattice  match.  In  the  case  of  quaternary  alloys,  there  are  two 
compositional  parameters,  so  there  is  an  additional  degree  of  free- 
dom which  can  be  exploited  to  obtain  an  optimum  band  gap  or  band 
lineup.  The  alloy  systems  which  are  currently  of  greatest  interest 
are  the  III-V  quaternary  systems  (Ga  ,In) (As , P)  and  (Al,Ga)(As,Sb) . 

Our  results  can  be  used  to  predict  the  band  lineups  of  hetero- 
junctions involving  alloys  by  using  a simple  interpolation  scheme. 
We  assume  that  the  band  edges  vary  linearly  for  the  ternary  alloys. 
This  does  not  quite  work  for  the  quaternary  alloys,  since  we  have 
four  pure  compounds  to  fit,  so  we  include  an  xy  term  to  account  for 
"warping"  of  the  band  edge  "surface."  For  the  GUj^Inj^.j^SyPj^.y 
alloy  system  this  scheme  gives  (in  eV) 

= -4.58  + 0.46x  + 0.20y  - 0.04xy 

(III.l) 


E^(r)  = -3.15+  1.89x  - 0.83y  - 0.35xy. 


The  corresponding  expressions  for  the  A1  Ga As  Sb,  system 

X A ” X y X y 

are 


E^  = -3.98  + 0.04x  + 0.02y  - 0.04xy 

E^(r)  = -3.04  + 1 .38x  + 0. 60y  - 0.  15xy  (III.  2) 

Ec(X)  = -2.81  + 0.46X  + 0.60y  - 0.25xy. 

This  phenomenological  interpolation  is  certainly  inaccurate 


to  the  extent  that  alloys  actually  exhibit  a small  nonlinearity  in  the 
band  gap  vs.  composition  curve.  We  have  attempted  to  calculate 
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the  alloy  band  structures  directly,  both  in  the  virtual  crystal 
approximation  and  with  a perturbative  correction  to  it  similar  to  that 
proposed  by  Baldereschi  and  Maschke  [1975] . However,  neither 
technique  successfully  reproduced  the  band  gaps  for  all  of  the  alloy 
systems  of  interest.  It  therefore  appears  that  it  will  require  a sig- 
nificantly different  approach  to  treat  the  alloys  properly,  a circum- 
stance already  noted  by  Van  Vechten  and  Bergstresser  [1970] . 

We  have  already  considered  the  GaAs-Al  Ga,  As  hetero- 

X 1 -X 

junction.  There  is  a great  deal  of  current  interest  in  the  InP- 
*^®0.48y^'^l-0.48y^^y^l-y  i'^^iction.  The  band  energies  as  a func- 
tion of  y may  be  evaluated  using  (Equations  (III.l).  The  results  are 
in  Figure  8.  Our  model  predicts  a very  small  conduction 
band  discontinuity  for  all  compositions.  In  some  very  preliminary 
experimental  results,  n-n  heterojunctions  of  this  system  have  been 
fabricated  and  they  show  nearly  ohmic  behavior  [R.  E.  Hayes, 
private  communication] , in  qualitative  agreement  with  our  predic- 
tion. 

The  predictions  for  heterojunctions  involving  the  III-V  quater- 
nary alloy  systems  (Ga  , In)(As  ,P)  and  (Al,Ga)(As  , Sb)  are  summarized 
in  Figure  9,  where  the  band  energies  are  plotted  as  a function  of 
lattice  constant.  Our  model  predicts  that  heterojunctions  within  the 
(Ga  , In)  (As  , P)  system  will  generally  have  most  of  the  discontinuity 
in  the  valence  band.  Those  within  the  (A1  ,Ga)  (As  , Sb)  system  will 
have  essentially  all  of  the  discontinuity  in  the  conduction  band. 


A 


(ln,Ga)As 


Figure  8.  Conduction  band  (upper  curve)  and  valence  band  (lower) 
curve)  energies  of  GaQ . 48y^'^l-0. 48y^^y^l-y  ^ 

function  of  composition  parameter  y. 


LATTICE  CONSTANT  (A) 


Figure  9 


. The  band  edge  energies  of  the  quaternary  alloy  systems 
(Ga  ,In) (As  ,P)  and  (A1  ,Ga)(As,Sb)  as  a function  of  lattice 
constant.  The  dotted  lines  indicate  indirect  conduction 
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Lattice-matched  (Ga  ,In)  (As  ,P) -(Al,Ga)  (As  ,Sb)  heterojunctions 

should  show  a stagger  in  band  lineup,  with  the  stagger  becoming  J 

especially  pronounced  for  alloys  with  large  lattice  constants. 

D.  Comparison  with  the  Electron-Affinity  Rule 

The  most  extensive  electron-affinity  data  available  for 
vacuum-cleaved  surfaces  are  those  of  Gobeli  and  Allen  [1966]  for 
Group  IV  and  III-V  semiconductors , and  those  of  Swank  [1967]  for 
the  II-VI  compounds.  It  is  important  to  use  only  data  for  vacuum- 
cleaved  surfaces,  because  the  electron  affinity  is  quite  sensitive 
to  any  surface  contamination. 

For  the  Ge-GaAs  heterojunction  the  electron-affinity  rule  pre- 
dicts a conduction  band  discontinuity  of  0.06  eV,  in  good  agreement 
with  both  our  prediction  and  the  experimental  data.  There  appear 
to  be  no  suitable  experimental  data  on  either  AlAs  or  GaP,  so  we 
cannot  compare  the  heterojunction  predictions  for  the  GaAs-ALAs  and 
Si-GaP  junctions.  The  electron  af  :inity  prediction  for  the  InP-CdS 
junction  (using  the  data  of  Fischer  [1966]  for  the  InP  electron 
affinity)  gives  a conduction  band  discontinuity  of  0.39  eV,  with  the 
InP  band  higher,  in  better  agreement  with  Shay,  et  al.  , [1976]  than 
with  our  prediction.  For  the  InP-(Ga , In)As  heterojunction  the  elec- 
tron affinity  rule  predicts  a conduction-band  discontinuity  of 
0.10  eV,  with  the  alloy  conduction  band  being  the  higher  one.  Our 


model  predicts  a discontinuity  of  the  same  magnitude,  but  in  the 
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opposite  direction.  Finally,  our  model  and  the  electron-affinity 
rule  are  in  qualitative  agreement  on  the  band  lineup  of  the 
InAs-GaSb  junction. 

It  therefore  appears  that  our  model  and  the  electron-affinity 
rule  generally  disagree  by  no  more  than  a few  tenths  of  an  electron 
volt.  This  level  of  disagreement  is  not  much  greater  than  the 
expected  accuracy  of  our  predictions  (and  the  expected  accuracy  of 
the  electron-affinity  rule  is  certainly  debatable)  so  it  is  not  clear 
that  the  differences  are  significant.  Our  model  has  the  advantage 
of  being  able  to  provide  a more  complete  set  of  predictions  than  the 
available  electron-affinity  data,  particularly  for  heterojunctions 
involving  GaP  and  AlAs.  It  is  also  more  readily  adapted  to  detailed 
studies  of  the  electronic  structure  of  the  interface. 


CHAPTER  IV 


DISCUSSION  AND  CONCLUSIONS 

Perhaps  the  most  significant  result  of  this  work  has  been  the 
demonstration  that  it  is  possible  to  define  a meaningful  "absolute" 
energy  scale  for  band  structure  calculations.  The  relationship 
between  the  band  structure  and  the  electrostatic  potential  is  cer- 
tainly unique,  if  somewhat  difficult  to  determine  computationally, 
and  the  success  of  our  simple  matching  scheme  demonstrates  the 
sort  of  understanding  that  can  be  gained  from  a knowledge  of  this 
relationship.  Our  calculations  are  certainly  not  the  last  word  on 
this  subject.  We  have  achieved  surprisingly  good  band  structures 
for  a model  with  only  two  adjustable  parameters  per  ion,  but  it 
would  be  very  desirable  to  apply  more  sophisticated  band  structure 
methods,  such  as  nonlocal  pseudopotential  and  OPW  methods,  to 
this  sort  of  calculation.  Perhaps  others  will  try  to  derive  the  elec- 
trostatic potential  from  their  calculations , now  that  we  have  shown 
that  this  is  useful  information. 

The  matching  scheme  deserves  a great  deal  more  attention 


and  development.  The  success  of  our  two-point  scheme  indicates 
that  the  physics  of  the  heterojunction  interface  is  likely  to  be 
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considerably  simpler  than  one  might  at  first  expect.  However,  it  | 

would  be  quite  desirable  to  have  a more  detailed  theory.  It  appears  j 

to  us  that  the  most  feasible  way  to  investigate  this  question,  ; 

j 

beyond  our  simple  scheme,  would  be  through  some  sort  of 
dielectric-response  model  of  the  detailed  charge  distribution  at  the 
interface.  Such  a scheme  would  presumably  be  sufficiently  simple 
that  it  could  be  applied  to  a wide  variety  of  heterojunctions.  On  a 
more  elaborate  level  one  could  apply  the  techniques  currently  being 
used  to  calculate  surfaces,  and  both  groups  presently  doing  such 
calculations  [Appelbaum,  1975;  Schluter,  1975]  have  expressed  an 
intention  to  do  just  that.  This  will  certainly  improve  our  under- 
standing of  heterojunction  physics,  but  it  will  not  invalidate  the 
simpler  models,  which  are  likely  to  contribute  more  to  our  under- 
standing of  such  aspects  as  chemical  trends  in  heterojunction 

i 

properties. 

With  our  present  two-point  matching  scheme,  our  model 
shares  the  property  of  transitivity  with  the  electron  affinity  model.  ; 

That  is,  if  we  know  the  band  lineups  for  heterojunctions  A-B  and 
B-C,  we  can  derive  the  lineup  for  A-C . Much  of  our  original 
objection  to  the  electron  affinity  rule  concerned  the  fact  that  there 
is  no  fundamental  reason  to  believe  that  heterojunctions  should 
display  this  property.  Only  a great  deal  of  further  experimental 


work  can  determine  the  extent  to  which  the  transitivity  property 
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actually  holds.  (Note  that  nontransitive  theories  necessarily 
involve  the  structure  of  the  interface  itself.) 

Our  results  can  also  be  viewed  as  part  of  a theory  of  electron 
affinities.  There  is  a correlation  between  our  conduction  band 
encgies  and  measured  electron  affinities,  but  there  is  a great  deal 
of  scatter  in  the  data  (Figure  10) . The  electron  affinity  should 
involve  three  contributions:  the  conduction  band  energy  with 
respect  to  the  electrostatic  potential,  the  surface  dipole  (definea 
as  the  difference  between  the  reference  potential  in  the  bulk  and  the 
vacuum  level)  , and  the  image  potential.  Our  calculations  give  the 
first  of  these  contributions.  As  our  understanding  of  surfaces 
improves,  it  may  be  possible  to  derive  a phenomenological  approxi- 
mation for  the  surface  dipole.  We  investigated  the  effect  of  image 
forces  at  heterojunctions  and  found  it  to  be  generally  negligible 
[Frensley,  1975],  but  this  is  probably  not  the  case  at  semiconductor 
surfaces . 

We  should  note,  at  this  point,  that  our  results  for  the  rela- 
tionship between  the  reference  potential  and  the  energy  bands 
depend  sensitively  on  the  choice  of  the  parameter  a in  the  local 
exchange  approximation,  although  the  relative  lineup  between  semi- 
conductors (which  is  what  is  important  for  heterojunction  lineups) 
is  much  less  sensitive.  Therefore,  any  attempt  to  relate  our  results 
to  absolute  quantities  such  as  the  electron  affinity  will  run  squarely 
into  the  issue  of  the  proper  choice  of  a.  Viewed  differently,  a 
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theory  of  electron  affinities  coulc  contribute  to  our  understanding  of 
the  appropriate  parameter  in  such  approximations,  though  it  would 
probably  not  contribute  to  a fundamental  understanding  of  the  local 
exchange  approximations. 

It  is  possible  that  our  calculations  could  contribute  to  the 
study  of  ionicity  of  semiconductors.  The  differences  in  the  electro- 
static potentials  at  our  reference  points  (AV) , as  given  in  Table  III, 
show  a good  correlation  with  ionicities  as  derived  from  electro- 
negativities (Figure  11).  The  correlation  with  the  Phillips  scale 
[Phillips,  1973,  p.  54]  is  somewhat  better  than  that  with  the 
Pauling  scale  [Pauling,  1960],  particularly  in  the  case  of  the  cad- 
mium compounds.  To  this  extent,  our  results  support  the  Phillips 
definition  of  electronegativity. 

There  is  another  approach  to  the  definition  of  a reference 
level  for  band  structures  which  deserves  mention.  McCaldin, 
McGill,  and  Mead  [1976]  have  used  the  Fermi  level  of  gold  in 
metal-semiconductor  junctions  as  a reference  level.  They  found 
that  the  valence  band  location  with  respect  to  this  level  is  only 
dependent  upon  the  anion  species  and  furthermore  correlates 
strongly  with  the  anion  electronegativity.  Our  results  show  a 
similar  correlation,  but  with  much  more  scatter  than  McCaldin, 
et  al.  repxirt.  This  subject  deserves  a great  deal  more  study,  par- 
ticularly to  determine  why  Schottky  barriers  using  copper  show  more 
scatter  than  those  using  gold,  and  how  this  is  related  to  the 
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detailed  electronic  structure  of  a metal-semiconductor  interface  as 
studied  by  Louie,  et  al.  [1976b]. 

Our  present  model  is  primarily  applicable  to  very  abrupt 
heterojunctions,  which  in  practical  terms  means  junctions  grown  by 
molecular  beam  epitaxy.  Liquid  phase  epitaxy  is  a bettor  developed 
technology,  which  produces  more  graded  heterojunctions.  The 
theory  of  graded  junctions  has  been  developed  by  Oldham  and 
Milnes  [1963]  and  a simplified  model  has  been  proposed  by  Cheung, 
Chiang,  and  Pearson  [1975].  These  models  stili  require  information 
about  the  "absolute"  band  energy  and  are  formulated  using  the  elec- 
tron affinity.  Our  results  for  the  conduction  band  energy  (assuming 
linear  interpolation  for  alloys)  could  be  used  in  these  models  in 
place  of  the  electron  affinity. 

In  conclusion,  we  have  developed  a model  for  the  band  lineup 
at  an  abrupt  semiconductor  heterojunction  which  is  remarkably  suc- 
cessful considering  its  simplicity.  The  heterojunction  lineup  is 
only  one  of  a wide  range  of  interface  phenomena  which  are  currently 
of  interest.  A semiconductor  hetero junction  is  one  of  the  "gentlest" 
possible  interfaces,  and  this  apparently  accounts  for  the  success 
of  our  model.  One  may  hope  that  the  understanding  we  have  gained 
can  be  applied  to  the  solution  of  other  interface  problems. 
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APPENDIX  A 


FREE  TON  SPECTRA 

The  ionic  pseudopotential  is  supposed  to  represent  the  effec- 
tive interaction  between  the  closed-shell  ion  and  an  electron.  ‘ 

Therefore,  a natural  test  to  apply  is  to  calculate  the  energy  levels  I 

of  an  electron  in  the  field  of  an  isolated  ionic  pseudopotential,  and 
compare  the  results  to  the  spectroscopically  observed  levels.  This 
test  is  very  useful  in  the  initial  choice  of  parameters,  and,  after  1 

we  have  fitted  the  band  structures,  a reasonable  fit  to  the  ionic  i 

i| 

levels  increases  our  confidence  that  we  have  chosen  a good  pseudo-  ; 

•| 

potential.  | 

Since  the  ionic  pseudopotential  which  we  are  using  is  a 

simple,  spherically  symmetric  potential,  calculation  of  the  elec-  j 

i 

1 

tronic  levels  is  a matter  of  solving  the  radial  part  of  the  j 

Schroedinger  equation.  This  is  conveniently  done  by  the  technique  ' 

of  finite-difference  equations.^  The  wavefunction  is  calculated  at 

a finite  number  of  points  and  the  differential  equation  is  transformed 

into  a tridiagonal  matrix  equation.  .j 

:i 

1 ! 

The  author  wishes  to  express  his  appreciation  to  i 

Mr.  David  G.  Unger  for  calling  this  technique  to  his  attention.  j 
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We  wish  to  find  the  eigenvalues  E of  the  radial  Schroeainger 
equation 

- 4-  f-  [ ^(r)  + V(r)v(r)  = E^(r)  (A.  1) 

r‘‘  dr  dr  y r 

where  V(r)  is  our  ionic  pseudopotential.  The  first  step  is  to  trans- 
form the  independent  variable  to 

X = . (A. 2) 

a+r 

Whereas  r runs  from  0 to  ■»,  x runs  from  0 to  1 . a is  effectively 
a scale  factor.  The  radial  equation  then  becomes 


- -Ti  + v[r(x)]4^(x)  = E^.(x) . 

3 X QX  cl  X 


(A.  3) 


If  we  define  f(x)=xH>(x)  , we  have  the  simpler  form 


f(x)  + V[r(x)]  = Ef(x)  (A.  4) 


dx 


Now  we  consider  a finite  mesh  of  n equally  spaced  points  x^=jA 
where  A=  l/(n+l) . The  f(x.]  form  an  n-dimensional  vector,  and  if 
we  approximate  the  derivative  by  finite  differences,  the  Hamiltonian 
becomes  a tridiagonal  matrix  of  order  n,  with  nonzero  elements 


H.  . = H.  , , 

J/J-1  J-J+1 


(1-Xj)^ 

a^A^ 


(A.  5) 
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This  matrix  can  then  be  diagonalized  to  find  the  eigenvalues  of  the 
original  Hamiltonian. 

A computer  program  was  written,  which  straightforwardly 
evaluated  the  expressions  (A.  5)  and  used  library  subroutines  to 
diagonalize  the  resulting  matrix.  A mesh  of  one  hundred  points  was 
used,  and  the  scale  factor  a was  set  to  2 a.u.  The  program  was 
tested  with  a Coulomb  potential  with  charge  Z=4 . The  calculated 
energies  were  within  0.003  Ry.  of  the  exact  values  for  all  states 
with  principal  quantum  number  n less  than  eight,  except  for  the  Is 
state,  which  receives  a significant  contribution  from  the  Coulomb 
singularity.  This  problem  will  not  occur  for  our  pseudopotential, 
since  it  is  smooth  and  continuous.  We  therefore  believe  that  the 
program  is  capable  of  accurately  calculating  the  energy  levels  of 
our  ionic  pseudopotential. 

The  calculated  ionic  spectra  for  the  pseudopotentials  actually 
used  in  the  band  structure  calculations  are  given  in  Table  IV,  along 
with  the  experimental  values.  The  calculated  values  generally 
agree  with  the  experimental  ones  to  within  about  5 percent,  which 
is  probably  about  as  good  as  one  could  expect  with  a local, 
energy- independent  pseudopotential.  The  values  for  the  Si  row  of 
the  periodic  table  illustrate  the  limits  of  a local  pseudopotential. 
The  d states  are  uniformly  lower  in  energy  than  the  calculated 
values.  This  is  because  there  are  no  occupied  d states  in  the 
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core,  so  the  d states  are  not  "pseudized,"  a point  discussed  by 
Chelikowsky  and  Cohen  [1974b]. 

The  importance  of  fitting  the  free  ion  spectra  is  illustrated  by 
the  fact  that  we  found  it  possible  to  obtain  reasonable  band  struc- 
tures for  the  antimonide  compounds  with  a different  pseudopotential 
than  that  which  was  finally  used.  This  potential  had  a narrower, 
taller  central  peak,  and  gave  a value  for  the  Sb  5s  state  which 
was  several  eV  too  low.  It  also  led  to  band  energies  with  respect 
to  the  reference  potential  which  were  about  0.3  eV  different  from  the 
final  values.  Therefore,  we  think  that  it  is  quite  important  to  apply 
such  a test  in  self-consistent  pseudopotential  calculations. 
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TABLE  IV 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDOPOTENTIAL 

The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


A1 


n 

s 

P 

d 

-28.45 

-21.77 

-14.07 

3 

-28.95 

-21.21 

-12.94 

-12.81 

-10.63 

-7.89 

4 

-13.27 

-10.46 

-7.31 

-7.29 

-6.32 

-5.03 

5 

-7.54 

-6.26 

-4.72 

-4.71 

-4.19 

-3.48 

6 

-4.86 

-4.17 

-3.30 

Si 

n 

s 

P 

d 

-45. 14 

-36.26 

-25.26 

3 

-46. 54 

-35.81  ■ 

-22.69 

-21.09 

-18.07 

-14.  14 

4 

-21.75 

-17.87 

-12.83 

-12.23 

-10.85 

-9.00 

5 

-12.57 

-10.76  - 

-8.29 

-7.99 

-7.24 

-6.22 

6 

-8.  19 

-7.21 

-5.82 
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TABLE  rv  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDO  POTENTIAL 

The  experimental  values  are  shown  directly  above- the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV  . 


P 

^ n 

s 

P 

d 

-65.02 

-53.97 

-39.71 

3 

-67.34 

-52.95 

-34.69 

-31.18 

-27.29 

-22.20 

4 

-32.31 

-26.85 

-19.69 

-18. 33 

-16.52 

-14.11 

5 

-18.91 

-16.32  - 

-12.76 

-12.06 

-11.08- 

-9.75 

6 

-12.41 

-10.99 

-8.96 

S T 

n 

s 

P * 

d 

.■*  , 

-88.05 

-74.82 

-57.37 

3 

-89.08 

-71.49'*_ 

-48.53 

-43.05 

-38.28 

-32.04 

4 

-44.03 

-36.97  _ 

-77  67 

-25.55 

-23.31«* 

-20. 35 

5 

-26.15 

-22.69 

--17.99  , 

-16.91 

-15.69 

-14.05 

6 

-17.30 

-15.37 

-12.67 
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TABLE  IV  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDO  POTENTIAL 

The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  ara  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


Zn 

n 

s 

P 

d 

-17.96 

-11.88 

-5.95 

4 

-17.79 

-11.37 

-5.99 

-7.00 

-5.38 

-3.34 

5 

-7.21 

-5.30 

-3.  37 

-3.77- 

-3.09 

-2.  14 

6 

-3.89 

-3.07 

-2.17 

-2.36 

-2.04 

-1.49 

7 

-2.43 

-2.01 

-1.52 

Ga 

n 

s 

P 

d 

-30.71 

-22.49 

-12.84 

4 

-32.42 

-22.67 

-13. 17 

-13.26 

-10.73 

-7.25 

5 

-14.29 

-10.97 

-7.43 

-7.46 

6 

-7.98 

-6.49 

-4.78 
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TABLE  IV'  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDO  POTENTIAL 


The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


Ge 


-45.71 

-35.40 

-22.06 

4 

-47.21 

-35.68 

-22.58 

-21.01 

-17.56 

-12.65 

5 

-22.25 

-17.88 

-12.79 

-12.23 

-10.54 

6 

-12.84 

-10.78 

-8.26 

As 

n 

s 

P 

d 

-62.63 

-50.24 

-33.22 

4 

-64.47 

-50.91 

-34.02 

-29.95 

5 

-31.64 

-26.19 

-19.38 

6 

-18.68 

-16.01 

-12.59 

i 
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TABLE  I\'  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDO  POTENTIAL 

The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


n 

s 

Se 

P 

i 

i 

\ 

d j 

4 

-81.7 

-83.98 

-67.3 

-67.84 

-46.6 
-47. 15 

5 

-40.3 

-42.75 

-35.74 

-27.06 

6 

-25.67 

-22.13 

-17.67  < 

n 

s 

Cd 

P 

d 

-16.91 

-11.23 

-5.78 

5 

-15.33 

-10.67 

-5.95 

-6.62 

-5.11 

-3.24 

6 

-6.51 

-5.06 

-3.35 

-3.60 

-2.97 

-2.09 

7 

-3.60 

-2.97 

-2.  16 

-2.27 

-1.95 

-1.46 

8 

-2.28 

-1,96 

-1.51 
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TABLE  IV  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDO  POTENTIAL 

The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


In 

n 

s 

P 

d 

5 

-28.03 

-28.02 

-20.58 

-21.18 

-12.08 

-12.99 

6 

-12.30 

-12.76 

-10.00 

-10.41 

-6.88 

-7.33 

7 

-7.03 

-7.29 

-5.60 

-6.22 

-4.46 
-4.  73 

8 

-4.55 

-4.73 

-4.15 

-3.31 

Sb 

n 

s 

P 

d 

5 

-55.7 

-57.11 

-44.9 

-46.31 

-30.6 

-32.44 

6 

-27.9 

-29.33 

-23.9 

-24.53 

-18.64 

7 


-17.69 


-15.24 


-12.20 


79 


TABLE  IV  (continued) 

FREE  ION  SPECTRA  DERIVED  FROM 
THE  IONIC  PSEUDOPOTENTIAL 

The  experimental  values  are  shown  directly  above  the  calculated 
ones.  Experimental  data  are  from  Moore,  1949,  1952,  and  1958. 
Energies  are  in  eV. 


' n 

s 

P 

d 

i 

[ 

* 

-72.3 

-59.8 

-42.7 

[ 5 

-74. 39 

-61.20 

-44.39 

1 

-37.8 

-33.0 

1 6 

-39.94 

-33.38 

-25.82 

7 


-24.55 


-21.03 


-17.01 


APPEJiDK  B 


COMPUTER  PROGRAMS 


1 . Band  Structure  Program 

The  primary  task  of  the  band  structure  program  is  to  calculate 
the  wavefunctions  and  energies  of  one-electron  states  in  a periodic 
potential.  We  therefore  make  use  of  Bloch's  theorem  [Kittel,  1963, 
ch.  9]  to  write  the  wavefunctions  as 

3 ° 


where  k is  a vector  from  the  first  Brillouin  zone,  the  S's  are 
reciprocal  lattice  vectors,  and  the  a-'s  are  complex  coefficients. 
For  computational  purposes  the  wavefunction  is  most  conveniently 
represented  by  the  one-dimensional  vector  of  the  sg's.  With  this 
representation  of  the  wavefunction  the  Hamiltonian  becomes  a 
matrix 


(B.2) 


The  first  term  is  the  kinetic  energy  (the  units  are  Rydbergs  and 
reciprocal  Bohr  radii,  which  are  used  for  all  internal  calculations). 
The  second  term  is  the  Fourier  coefficient  of  the  potential  for  the 
reciprocal  lattice  vector  S-5'. 


81 


At  this  point  it  should  be  obvious  that  in  order  to  represent 

our  states  and  Hamiltonian  as  subscripted  arrays,  we  must  number 

the  reciprocal  lattice  vectors  so  that  they  may  be  used  as  indicies. 

This  is  done  by  the  subroutine  FCCLAT,  which  constructs  the 

reciprocal  lattice  vectors  for  the  1 sce-centered  cubic  lattice.  All 
2 

vectors  with  G = 32  are  constructed,  sorted  in  order  of  increasing 
2 

S , and  stored  in  the  TABLE  common  block.  A cross-table  which 
will  give  an  index  when  entered  with  the  three  components  of  S is 
then  constructed  and  stored  in  TABLE. 

The  potential  is  constructed  by  either  FRMFCTR  or  MODELV, 
depending  upon  whether  an  empirical  pseudopotential,  such  as 
those  of  Cohe.n  and  Bergstresser  [1966],  or  the  self-consistent 
potential  described  in  Chapter  II  is  to  be  used.  It  is  convenient  to 
choose  the  origin  of  the  coordinate^ system  at  the  point  halfway 
between  the  two  atoms  in  the  unit  cell.  The  potential  can  then  be 
written  as 

v(G)  = V (G)cosG‘T  + iv  (G)  sinG'T  (B.3) 

s a 

where  t =a  , Vg  and  v^  are  the  symmetric  and  antisymmet- 

ric parts  of  the  potential  under  interchange  of  the  two  atoms. 
Vg{G)=0  for  the  group  IV  semiconductors. 

The  actual  construction  and  fliagonalization  of  the 
Hamiltonian  is  done  by  the  subroutine  EVAL.  This  routine  is  entered 
with  the  k vector  of  the  desired  Bloch  states . EVAL  calculates  the 


kinetic  energy  of  the  basis  states  Ik+G)  and  calls  the  SORT  routine 
which  sorts  the  states  in  order  of  increasing  energy.  The  basis 


states  are  then  partitioned  into  the  two  sets  required  by  the  Lowdin 

2 

procedure.  The  first  set  includes  those  states  with  I k+G  I =E^ 

2 

and  the  second  includes  those  with  Ej^<  I k+G  I =^2’  then 

constructs  the  parts  of  the  Hamiltonian  required  by  the  Lowdin 
procedure,  and  calls  LOWDIN3,  which  iiarforms  that  procedure,  as 
mc'uified  by  Brust  [1964].  The  higher-energy  states  are  "folded 
down"  as 


H 

mn 


t E-Hjj 


(B.4) 


where  m and  n are  states  from  the  first  set,  j is  a state  from  the 
second  set,  and  E is  either  the  average  energy  of  the  first  eight 
states,  for  m/n,  or  H for  n=m.  EV^L  then  calls  library  matrix 
diagonalization  routines . __ 

If  the  eigenvectors  are  required,  caller*  set  of  routines 

which. perform  that  calculation.  It  then  restores  the  folded-down 
components  by  . 

The  components  must  finally  be  rearranged  so  as  to  conform  to  the 
original  numbering  of  the  states  as  given  in  TABLE. 

The  calculation  of  the  valence  charge  distribution  is  super- 

« 

vised  by  the  subroutine  VALCiIG.  It  sets  up  the  k vectors  required 
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for  the  Chadi-Cohen  [1973a]  three-point  Bnllouin  zone  sample  and 
then  calls  EVAL  to  calculate  the  eigenvectors.  The  subroutine 
ABSPSI2  calculates  the  absolute  square  of  the  Bloch  function  by 


! v(x)  I ^ 


G,G' 


i(G-G')-x 


16.6) 


and  then  adds  this  to  the  charge  density  array  in  common  block 
CHGDEN,  with  the  appropriate  multiplicative  factor.  After  all  of 
the  wavefunctions  have  been  calculated  and  squared,  VALCHG  calls 
SYMETRZ  to  symmetrize  the  charge  density  with  respect  to  the  point 
group  operations.  Since  we  have  chosen  the  origin  of  our  coordi- 
nate system  at  a point  which  does  not  exhibit  the  full  point  symme- 
try of  the  crystal,  we  must  first  translate  the  origin  to  one  of  the 
atomic  sites.  This  is  done  by  multiplying  the  Fourier  coefficients 
of  the  charge  density  by  a phase  factor.  The  subroutine  then  runs 
through  the  24-point  group  operations,  which  are  represented  as 
3X3  matricies  R,  operating  on  each  G and  adding  the  G coeffi- 
cient of  the  charge  density  to  the  RG  coefficient  of  the  symme- 
trized charge  density  ar.ay.  The  coordinate  system  must  be  finally 
translated  back  to  the  midpoint  of  the  bond,  again  by  multiplication 
by  phase  factors.  If  desired,  the  Fourier  series  for  the  charge 
density  can  be  summed,  and  the  results  plotted  out  along  the  body 


diagonal,  by  the  subroutine  FSPLT. 
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The  Slater  approximation  for  the  exchange  interaction  is 
evaluated  from  the  charge  density  by  SLATERS.  The  quadratic 
approximation  (II. 9)  for  the  cube  root  is  used.  The  square  of  the 
charge  density  is  calculated  by  CONVOL,  which  performs  the  dis- 
crete convolution  needed  to  find  the  product  of  two  Fourier  series , 
given  by 

y , . iS’x  V K-  ^i(5+S’)*x  , . 

2,(ab)ge  = I 

S G,G' 

There  are  three  subroutines  which  supervise  the  execution  of 
the  major  tasks  the  program  is  intended  to  perform.  ITERATE  per- 
forms the  self-consistency  iterations.  Using  the  current  valence 
charge  distribution,  ITERATE  calls  MODELV  to  calculate  the  pseudo- 
potential. It  next  calls  VALCHG  to  perform  the  eigenvector  and 
charge  density  calculations.  The  resulting  charge  density  is  taken 
in  linear  combination  with  the  previous  density,  in  order  to  provide 
a means  for  controlling  oscillatory  convergence  of  the  band  struc- 
ture. Finally  the  average  root-mean-square  change  in  the  electro- 
static potential  is  calculated  and  printed  out.  The  subroutine 
PRINGAP  is  used  to  calculate  the  principal  gaps  at  the  points  r, 

L,  and  X.  It  calls  EVAL  to  calculate  the  appropriate  eigenvalues, 
assigns  symmetry  labels  to  the  levels  at  r on  the  basis  of  the 
degeneracies,  and  prints  out  the  resulting  gaps.  BZSAMPL  does  a 
more  complete  sampling  of  the  Brillouin  zone,  calculating  the  bands 


along  the  A,  A,  and  Z directions.  It  calls  the  subroutine 
BANDPLT,  which  is  a straightforward  printer  plotting  routine,  to 


plot  the  results.  i 

The  main  program  BANDST  performs  none  of  the  actual  calcu- 
lations. It  handles  the  disk  operations,  modifies  the  parameter  set 
in  response  to  input  data,  decodes  input  instructions,  and  calls  the  i 

subroutines  to  perform  the  calculations.  It  uses  a number  of  fea- 
tures peculiar  to  CDC  systems  and  tne  University  of  Colorado  com- 
puting facilities  in  particular,  so  it  would  certainly  have  to  be 

i 

extensively  rewritten  in  order  to  run  successfully  elsewhere. 

2.  Electrostatic  Potential  Program 

The  electrostatic  potential  is  calculated  according  to  the 
scheme  described  in  Chapter  II.  This  involves  subtracting  a ficti- 
tious charge  distribution  from  the  valence  charge  and  adding  it  to 
the  ion  cores.  The  subroutine  FOURCF  subtracts  the  fictitious  dis- 
tribution (II. 8)  from  the  valence  charge  and  divides  the  Fourier 
coefficients  by  to  find  the  electrostatic  potential.  The  Fourier 
series  for  this  part  of  the  potential  may  be  summed  either  by  FSPT, 
if  the  value  at  only  one  point  is  desired,  or  by  FRSUM , if  the 
values  along  a line  through  the  crystal  are  desired. 

The  ionic  contribution  to  the  electrostatic  potential  is  cal- 
culated by  the  subroutine  CLUSTER.  It  constructs  the  lattice 
points  of  the  face-centered  cubic  structure,  calculates  the 


F 
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distances  from  the  given  point  to  the  two  atomic  sites,  and  calls 
the  appropriate  atomic  potential  function,  VATOMA,  etc.  with 
these  distances.  It  then  adds  the  resulting  contributions  onto  the 
total  potential. 

3.  Program  Listings 

These  programs  have  been  written  for,  and  successfully  run 
on,  the  CDC  6400  system  at  the  University  of  Colorado,  operating 
under  KRONOS  2. 1 and  using  the  RUN  Fortran  compiler. 

As  listed  in  the  following  pages,  the  program  CRYSPOT  and 
its  associated  subroutines  assume  locations  for  the  anion  and  cation 
within  the  unit  cell  which  are  reversed  from  those  locations  assumed 
by  BANDST  and  its  associated  routines.  The  following  changes  will 
make  the  two  programs  consistent.  Note  that  this  affects  only  the 
interchange  of  information  between  the  programs;  each  will  run  suc- 
cessfully as  listed. 

Subroutine  FOURCF  Line  above  statement  no.  200  should  be 
AR(I)  =-{C/GSQ)*(CI(D  - (RHOA-RHOB)*SINGTAU(D) 

Subroutine  CLUSTER 

Statement  no.  100  should  be 

100  R2  =R2+ (Y(I) -TAU(I) -X(D)**2 

Statement  no.  170  should  be 

170  R2  =R2+ (Y(I) +TAU(I)  -X(I))**2 
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PROGRAM  BANDST  (INPUT.OUTPUT, PUNCH. TAPE1*1NPUT. TAPE?) 


LOGICAL  INFOBIT 


C0ft10N/INF0/INF0BIT(16),El,E2,A 

THE  INFO  COMTTON  BLOCK  CONTAINS  LOGICAL  VARIABLES  CONTROLLING 

THE  OPERATION  OF  THE  PROGRAM  AND  THE  OUTPUT 

INFOBIT(I)  PRINTS  OUT  TIMING  INFORMATION 

INF0BIT(2)  PRINTS  OUT  MILLER  INDEX  TABLE 

INF0BIT(3)  PRINTS  OUT  EIGENVALUES 

INF0BIT(4)  PRINTS  OUT  CHARGE  COMPONENTS 

INF0BIT(5)  PRINTS  OUT  CHARGE  PLOT 

INF0BIT(6)  INDICATES  THAT  AN  EPH  PSEUOOPOTENTIAL  IS  TO  BE  USED 
El  IS  THE  LIMIT  OF  THE  INNER  LOWDIN  SET 
E2  IS  THE  LIMIT  OF  THE  OUTER  LOUDIN  SET 
A IS  THE  LATTICE  CONSTANT 


COMMON/LOCALV/VR( 181 ) . VI ( 1 81 ) 

THE  LOCALV  COMflON  BLOCK  CONTAINS  THE  FOURIER  COMPONENTS  OF  THE 
LOCAL  POTENTIAL. 


C0»WN//VMNR(I600)  ,VMNI(1600)  ,VJMR{3200)  ,VJMI  (3200)  ,VJJ(80) 

THE  BLANK  COMMON  BLOCK  CONTAINS  THE  RELEVANT  PARTS  OF  THE 
UNPERTURBED  HAMILTONIAN.  THE  VMN  CONTAIN  THE  SMALL  SQUARE 
MATRIX.  THE  VJM  CONTAIN  THE  RECTANGULAR  MATRIX.  AND  VJJ 
CONTAINS  THE  DIAGONAL  ELEMENTS  OF  THE  LARGE  SQUARE  MATRIX. 


COMMON/HAMLT/HR ( 1 600 ) .H I ( 1 600 ) 

THE  HAMLT  C0MI40N  BLOCK  CONTAINS  THE  REAL  AND  IMAGINARY  PARTS  OF 
THE  PERTURBED  HAMILTONIAN. 


COMMON/TABLE/ ITAB( 3 . 1 81 ) .NTAB( 1 1 . 1 1 , 1 1 ) 

THE  TABLE  COMMON  BLOCK  CONTAINS  THE  INFORMATION  REQUIRED  TO 
CONVERT  FROM  MILLER  INDICIES  TO  A SINGLE  SUBSCRIPT  AND 
VICE  VERSA. 


COWON/HEAD I N6/  NAME  ( 8 ) . I DATA  ( 8 ) 

THE  HEADING  COMMON  BLOCK  CONTAINS  CHARACTERS  USED  IN  OUTPUT 
HEADINGS. 


COMMON/PTGROUP/GKN (3.3.24) 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 


THE  PTGROUP  COWlON  BLOCK  CONTAINS  THE  3X3  flATRICIES  OF  THE 
OPERATIONS  OF  THE  POINT  SYMMETRY  GROUP. 


COMMON/SFACTR/COSGTAU( 181 ) ,SINGTAU( 181 ) 

THE  SFACTR  COMMON  BLOCK  CONTAINS  THE  STRUCTURE  FACTORS 
COS(G*TAU)  AND  SIN(G*TAU)  FOR  THE  ZINC-BLENDE  STRUCTURE. 

C0MM0N/BAN0S/EK(8,40) 

THE  BANDS  COMMON  BLOCK  CONTAINS  THE  ENERGIES  (IN  EV)  AT  EACH  K. 
COMMON/CHGOEN/CR{ 181 ) ,CI{181 ) 

THE  CHGDEN  COMMON  BLOCK  CONTAINS  THE  FOURIER  COEFFICIENTS  OF 
THE  VALENCE  ELECTRON  CHARGE  DISTRIBUTION. 


C0MM0N/INCLUDE/INC(181) 

THE  INCLUDE  COMMON  BLOCK  CONTAINS  INFORMATION  ^COMPLEMENTED) 
AS  TO  WHICH  PLANE  WAVES  WERE  INCLUDED  IN  THE  EIGENVECTOR 
CALCULATION. 


C0M40N/PARAMTR/ZA,QA,ALPHAA,BETAA.GAMMAA,V0A,ZB,QB.ALPHAB,BETAB, 
IGAMMAB.VOB.XALPHA, RELAX 

C0MM0N/PRMTR2/ZC,QC,ALPHAC,BETAC,GAMMAC,V0C,ZD.QD,ALPHAD,BETAD, 

IGAfWAD.VOD.Xl.YI 

THE  PARAMTR  COLtlON  BLOCK  CONTAINS  THE  PARAMETERS  OF  THE  MODEL 
POTENTIAL. 


C0ft10N/EPMFF/NFF,IG?(10),VS(10),VA(10) 

THE  EPMFF  COItlON  BLOCK  CONTAINS  THE  PSEUDOPOTENTIAL  FORM 
FACTORS  FOR  USE  IN  EPM  CALCULATIONS. 


DIMENSION  IMAGE ( 80 ) . I NSTRCT ( 1 0 ) , I TEMP (2 ) 
DIMENSION  PRMTR(14) 

DIMENSION  PRMT2(14) 

EQUIVALENCE  (PRMTR{ 1 ) ,ZA) 

EQUIVALENCE  (PRMT2(1 ) ,ZC) 


90 


I 


CALL  GET  (5HTAPC7,7HPTGROUP,0,0) 

READ  (7,50)  {Gf1N(I).I  = 1.216) 

50  FORMAT  (9F5.0) 

REWIND  7 
C 
C 

CALL  GET  (5HTAPE7,7HPARAMTR.0,0) 

READ  (7,60)  (PRMTR(I),I=1,14),E1,E2,A 
READ  (7,61)  (NAME{1),I*1,8),(1DATA(I),I=1,8) 

READ  (7,60)  (PRHT2(I),I=1,14) 

60  FORMAT  (F12.4) 

61  FORMAT  (8A10) 

REWIND  7 

C 

C 

C READ  AND  ENCODE  COMMANDS 

C 

READ  100,(IMAGE(1),I=1 ,80) 

100  FORMAT  (80R1) 

C 

IP0INT=0 

C 

DO  120  J=1,10 
C 

I1=IP0INT+1 

C 

C 

DO  no  1=1,11 
c 

IP0INT=IP0INT+1 

C 

IF  (IMAGE(IPOINT).EQ.IR,)  GO  TO  115 
IF  (IMAGE(IPOINT).EQ.IR.)  GO  TO  125 
C 

110  CONTINUE 
C 

115  NCHAR=IP0INT-I1 
IENO=IPOINT-1 

ENCODE  (UCHAR,118,INSTRCT(J))  (IMAGE(1A),IA=I1  ,IEND) 
118  FORMAT  (lORl) 

C 

C 

120  CONTINUE 
C 

125  NCHAR=IP0INT-I1 
IEND=IP0INT-1 

ENCODE  (NCHAR,118.INSTRCT(J))  (IMAGE(IA),IA=I1 ,IEND) 
J*J+1 

INSTRCT(J)=10HHALT 

C 

C 

C READ  AND  ENCODE  OUTPUT  INSTRUCTIONS 

C 


c 


READ  100, ( IMAGE (1), I =1,60) 
DO  130  1=1,16 


130  INrOBlT(l)-.F. 

C 

11  = 1 
C 

00  140  1=1,80 
C 

IF  ((lMAGE(l).fiE.lR,).ANO.(inAGE{I).N£.lR.))  GO  TO  140 
C 

NCHAR=1-I1 

IEN0=I-1 

ENCODE  (NCHAR.llS.IWORD)  ( IMAGE( lA) . IA=I 1 . lENO) 

11=1+1 

C 

IF  (IWORD.FQ.lOHTIfllNG  ) INFOBIT(  1 ) = . T. 

IF  IWORD.FQ.IOHTABLE  ) INF0BIT(2)=.T. 

IF  (iWORD.EQ.lOrlEMERGIES  ) 1NF0B1T{3)  = .T. 

IF  (IWORD.EQ.lOhCHARGECOEF)  IGF03:T(4)=.T. 

IF  (IWORD.EQ.IOHCHARGEPLOT)  INF0BIT(5)=.T. 

C 

IF  (IMAGE(I).EQ.in.)  GO  TO  145 
C 

140  CONTINUE 
C 
C 
C 

C READ  HEADING, DATA, AND  FORM  FACTORS. 

C 

145  READ  150,IW0RD 
C 

150  FORMAT  (6A10) 

C 

IF  (E0F,1)  170,155 
C 

155  IF  (IWORD.EQ.IOHHEADING  ) READ  150, (NAME{ I ) , 1=1 ,8) 

IF  (IWORO.EQ.IOHOATA  ) READ  1S0,(I0ATA(1) ,1=1 ,8) 

C 

IF  (IWORD.NE.IOHFORM  FACTO)  GO  TO  165 
C 

READ  160, NFF 

READ  160,{IG2(I),VS(1),VA(I),I=1 ,NFF) 

150  FORMAT  (I5,2F10.4) 

C 

165  GO  TO  145 
C 
C 

C READ  MODEL  PARAMETERS 

C 

170  READ  100,(1MAGE(I),I=1,80) 

C 

IF  (E0F,1)  200,175 
C 

175  DO  180  1A=1 ,11 
C 

IF  (IMAGE(IA).EQ.1R=)  GO  TO  185 
C 

180  CONTINUE 


oooooo  oooo 


c 

c 

185  NCHAR=IA-1 
C 

ENCODE  (NCHAR,118,IW0R0)  {UIAGE(I)  ,1  = 1 .NCHAR) 
C 

IST=IA+1 
1END=IST+19 

ENCODE  (20,118.ITEMP)  ( IMAGE (’). I=1ST, lEND) 
DECODE  (20,190, ITEMP)  X 
190  F0RI-!AT  (F20.6) 


TEST  NAMES  AND  REPLACE  DATA 


IF 

(IW0RD.EQ.10H2A 

) ZA=X 

IF 

nwORO.EQ.lOHZB 

) ZB=X 

IF 

(IWORD.EQ.IOHQA 

) QA=X 

IF 

IWORD.EQ. lOHOB 

QB=X 

IF 

(IWORO.EQ. lOHALPHAA 

) ALPHAA=X 

IF 

(IWORD.EQ. lOHALPHAB 

) ALPHA6=X 

IF 

(IWORD.EQ. lOUBETAA 

) B£TA/>=X 

IF 

(IWORD.EQ.  lOtlBETAB 

) BETAB=X 

IF 

(IWORD.EQ.IOHGAMMAA 

) GAMMAA=X 

IF 

(IWORD.EQ.IOHGAMMAB 

) GAMMAB=X 

IF 

(IWORD.EQ. lOHVOA 

) VOA=X 

IF 

(IWORD.EQ. lOHVOB 

) VOB=X 

IF 

(IWORD.EQ. lOHA 

) A=X/.529177 

IF 

(IWORD.EQ. lOHEl 

) E1=X 

IF 

(IWORD.EQ. 10HE2 

) E2=X 

IF 

(IWORD.EQ. lOHXALPHA 

) XALPHA=X 

IF 

(IWORD.EQ. lOHRELAX 

) RELAX=X 

IF 

(IWORD.EQ. lOHZC 

) ZC=X 

IF 

(IWORO.EQ. lOHZD 

) ZD=X 

IF 

(IWORD.EQ.IOHQC 

) QC=X 

IF 

(IWORD.EQ.IOHQD 

) QD=X 

IF 

(IWORD.EQ. lOHALPHAC 

) ALPHAC=X 

IF 

(IWORO.EQ. lOHALPHAO 

) ALPHAD=X 

IF 

(IWORD.EQ. lOHBETAC 

) BETAC=X 

IF 

(IWORD.EQ. lOHBETAD 

) BETAD=X 

IF 

(IWORD.EQ. 10HGAMMAC 

) GAMMAC=X 

IF 

IWORD.EQ. lOHGAMMAO 

) GAMMAD=X 

IF 

(IWORD.EQ. lOHVOC 

VOC=X 

IF 

(IWORD.EQ. lOHVOD 

) VOD=X 

IF 

(IWORO.EQ. lOHX 

) X1=X 

IF 

(IWORD.EQ. lOHY 

) Y1=X 

GO 

TO  170 

PRINT  HEADING 

200  CALL  DATE  (II) 
CALL  TYME  (12) 


r 


9i 


c 

c 

c 

c 


c 

c 


c 


c 

c 

c 


c 

c 


c 

c 

c 

c 


PRINT  ?ns,!l.I2 

PRINT  ?’.0,(NAI!i;(I).I  -1  ,tt) 

205  FORMAT  (•  RUN  DATE  * . AlO , lOX . *T !«f  *,A10,/) 
210  FORMAT  (•  UAND  STRUCTURE  FOR  *,8A10,//) 


PRINT  POTENTIAL  PARAMETERS 

PRINT  215.F1.E2 
PRINT  220, XALPhA, RELAX 

A1=.529177*A 


PRINT  225, A1 
PRINT  230 

PRINT  235.ZA,7B.ZC,Z0 

PRINT  236.QA.C0.QC.0n 

PRINT  237  . Al  PI.AA.ALPHAB.ALPHAC.ALPHAD 

PRINT  238, BETAA. BETAS. BLTAC.5ETAD 

PRINT  239 . r.;.T1A A . OA.MVAO . GAMMAC . GAMMAD 

PRINT  2A.0.V0A.V0B.V0C.V00 

PRINT  245.X1.Y1 


215  FORMAT 
220  FORT’AT 
225  FORMAT 
230  FORMAT 

235  FORMAT 

236  FORMAT 

237  FORMAT 

238  FORI'AT 

239  FORMAT 

240  FORMAT 
245  F0RI4AT 


(*  CALCULATIONAL  PARAMETERS*. ET =* .F3.0 . 1 1 X .*E2=* .F3.0) 
(»  XALPHA=*.F4. 2. 5X, ‘RELAXATION  FACT0R=* .F4 . 2) 

(//.*  POTENTIAL  PARAMETERS*.//.*  LATTICE  CONSTANT=* ,FG. 3 ./ ) 
(TAX. ‘ATOM  A*,9X,*AT0M  B*,9X,*AT0M  C*.9X.*AT0M  D*./) 
(9X.*Z*.4(F10.0.5X)) 

(9X.*Q*.4(F10.0.5X)) 

(5X.*ALPHA*,4(F10.2.5X)) 

{6X.*BETA*.4(F10.2.5X)) 

(5X.*GAMMA*.4(F10.2.5X)) 

(8X.*V0*.4(FI0,2,5X)) 

(//,8X,*X=*,F4.2,6X,*Y=*,F4.2,///) 


WRITE  (7,60)  (PPMTR(I),I  = 1 ,14),E1  ,E2,A 
WRITE  (7,61)  {NAME(I),I=T,8),(IDATA(I),I=T,8) 
WRITE  (7,60)  (PRMT2(I),I=T  ,14) 

CALL  REPLACE  (5HTAPE7 ,7UPARAMTR,0,0) 

REWIND  7 


CALL  GET  (5MTAPE7,7HCMGDATA,0,0) 
READ  (7.20G)  (CR( I ) ,CI ( I ) , 1 = 1 ,T8T ) 
295  FORMAT  (2F10.5) 

REWIND  7 


EXECUTE  THE  TASKS  SPECIFIED  IN  INSTRCT 


300 


CALL  FCCLAT 


C 

c 


DO  450  1NST=1 ,10 


c 

C IMPLfMfNT  THE  INITIALIZE  INSTRUCTION 

C 

IF  (INSTRCT(INST).NE.IOHINITIALIZE)  go  to  350 
C 

CALL  FRflFCTR 
CALL  VALCHG 


C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 


C 

C 

C 


C 


C 

C 

C 


C 

C 

C 

c 

c 


irPLEMENT  THE  EPM  INSTRUCTION 

350  IF  (INSTRCT(INST).EQ.IOHEP"  ) INF0BIT{6)=.T. 

IMPLEMENT  THE  GAPS  INSIRUCTION 

IF  (INSTRCT(IMST).NE.10HG/'.rS  ) GO  TO  370 

IF  (INF0BIT(6))  CALL  FRMFCTR 
IF  (.N0T.INF0BIT(6))  CALL  MOOELV  (CR.CI) 

CALL  PRINGAP 

implement  the  full  zone  INSTRUCTION 

370  IF  (INSTRCT(INST).NE.IOHFULL  ZONE  ) GO  TO  390 

IF  (INF0BIT(6))  CALL  FRMFCTR 
IF  {,N0T.1NF0B1T(6))  CALL  MODELV  (CR.CI) 

CALL  BZSAMPL 

IMPLEMENT  THE  ITERATE  N INSTRUCTION 

390  I1=INSTRCT(INST).AND.077777777777777000000 
IF  (I1.NE.7LITERATE)  GO  TO  420 

I2=( INSTRCT{ INST) . AND. 0777777) .OR. 7L 
DECODE  (9,400,I2)NI 
400  FORMAT  (7X.I2) 

CALL  ITERATE  (Nl.RELAX) 

IMPLEMENT  THE  HALT  INSTRUCTION 

420  IF  (INSTRCT(INST).EQ.IOHHALT  ) GO  TO  500 

450  CONTINUE 

500  CONTINUE 


STORE  CHARGE  DENSITY  ON  DISK 

WRITE  (7,295)  (CR( I ) ,CI ( I ) , 1=1 . 181 ) 
CALL  REPLACE  (5HTAPE7 ,7HCHG0ATA,0,0) 


C 

C 

c 


END 
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SUBROUTINE  RCLAT 
C 
C 
C 

C THIS  SUBHOUTINE  SETS  UP  THE  TABLE  OF  RECIPROCAL  LATTICE  VECTORS 

C FOR  THE  FCC  STRUCTURE  AND  CALCULATES  THE  STRUCTURE  FACTORS. 

C 

C 

C0f»10N/TABLE/ITAB{3,181  ),NTAB(11  ,11,11) 

COMMON/SFACTR/COSGTAU( 181 ) ,SINGTAU( 181 ) 

C0FLM0N/INF0/INF0BIT(16) 

C 

c 

DIMENSION  1GSQ{181) 

C 

INTEGER  GSQ 
C 
C 

C SET  UP  TABLE  TO  CONVERT  MILLER  INDICIES  TO  SINGLE  SUBSCRIPT 
C 

C MA=1  CALCULATES  ODD  MILLER  INDICIES,  f1A=2  CALCULATES  EVEN  ONES 

C 

N = 1 

DO  210  MA=1 ,2 
MAX=5 

M1N=MA-MAX-1 

C 

C THESE  LOOPS  CONSTRUCT  THE  G VECTORS 

C 

DO  210  I=MIN,MAX,2 
00  210  J=MIN,MAX,2 
DO  210  K=MIN,MAX,2 
C 

GSQ=I*I+J*J+K*K 
IF  (GSQ.GT.32)  GO  TO  210 
C 

C THIS  SORT  ROUTINE  ORDERS  THE  VECTORS  IN  INCREASING  MAGNITUDE 

C 

L=N-1 

IF(N.EQ.l)  GO  TO  206 

204  IF(L.LE.O)  GO  TO  206 
IF(IGSQ(L).LE.GSQ)  GO  TO  206 
IGSQ(L+1)=IGSQ(L) 

DO  205  IA=1 ,3 

205  ITAB(IA,L+1)=ITAB{IA,L) 

L=L-1 

GO  TO  204 

206  IGSQ(L+1)=GSQ 
ITAB(1  ,L+1)  = I 
ITAB(2,L+1)=J 
ITAB{3,L+1)=K 
N=N+1 

210  CONTINUE 
C 

C SET  UP  TABLE  OF  N AS  A FUNCTION  OF  MILI  ER  INDICIES 

C 
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DO  212  N=1 .181 
IMTABd  ,\) 

J=ITAfi(2,fO 

K=ITAB(3,N) 

NTAB(W6,J+6,K+6)  = N 
212  CONTINUE 
C 
C 

C PRINT  OUT  INDEX  TABLE  IF  INFOBIT(2)  IS  TRUE. 

C 

IF  (.NOT.INFOBIT(2))  NO  TO  280 
C 

PRINT  250 
PRINT  255 

PRINT  260,(I.(1TAB(J,I),J=1.3),IGSQ(I),;=1,181) 

C 

250  FORMAT  (/// ,42X ,*TABLE  OF  RECIPROCAL  LATTICE  VECTORS*C//) 
255  FORMAT  (5{*  INDX  ( I J K)  GSQ*.7X)./) 

260  FORMAT  (5(15,*  (*,312,*)  *,I3,7X)) 

C 

C 

C SET  UP  STRUCTURE  FACTOR  TABLE. 

C 

C 

280  DO  300  1=1 ,181 
C 

GTAU=0.785398*(ITAB(1 ,1)  + ITAB(2,I)  + ITAB(3,D) 

C 

C0SGTAU(I)=C0S(GTAU) 

S1NGTAU(1)=SIN(GTAU) 

C 

C 


SUBROUTINE  FRMFCTR 
C 

C THIS  SUBROUTINE  CALCULATES  THE  FOURIER  COMPONETS  OF  THE 

C PSEUDOPOTENTIAL  FROM  THE  FORM  FACTOR  INFORMATION  STORED  IN  E 

C EPMFF. 

C 

r 

C0MM0N/EPMFF/NFF.1G2(10),VS(10).VA(10) 

COMMON/SFACTR/COSGTAU( 181 ) ,SINGTAU( 181 ) 

COMMON/LOCALV/VR( 1 81 ) . VI ( 1 81 ) 

COMMON/TABLE/IT{3, 181 ),NT( 11,11 .11) 

C 

C 

DO  100  1 = 1,181 
VR(I)=0. 

100  VI(I)=0. 

C 

C 

DO  200  1=2,181 
C 

IGS0=O 

DO  120  J=l,3 

120  IGSQ=1GSQ+IT(J.1)**2 
C 

DO  150  J=1,NFF 
C 

IF  (IGSQ.NE,IG2(J))  GO  TO  150 
C 

VR(I)=C0SGTAU(I)*VS(J) 

VI(I)=SINGTAU(I)*VA(J) 

C 

150  CONTINUE 
200  CONTINUE 
C 

RETURN 

END 


o o o <->  o r>  o o o oo  ooo  t~>  o ooooo 
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SUBROUTINE  MODELV  (CR.CI) 

THIS  SUBROUTINE  CALCULATES  THE  POTENTIAL  FROM  THE  MODEL 
PARAMETERS, 


C0W0N/TABLE/IT(3.181),NT(1 1.1 1,11) 

COIt10N/L0CALV/VR(  181 ) ,VI  ( 1 81 ) 
C0W^0N/SFACTR/C0SGTAU(181),SINGTAU(181) 

COMMON/ INFO/DUMMY (18), A 

COPWON/PARAfITR/  ZA  ,QA , AL  PHAA  ,8ETAA , GAMMAA , VOA , ZB  ,QB , AL  PHAB , BETAB , 
1 GAMMAB , VOB , XAL  PHA , RE  L AX 

C0MM0N/PRMTR2/ZC,QC,ALPHAC,BETAC,GAMMAC,V0C,ZD,QD,ALPHAD.BFTAD, 

1GAMMA0,V00,X1,Y1 

DIMENSION  CR(181),CI{181) 

REAL  KSQ 


CALL  SLATER3  (XALPHA, CR.CI ,VR, VI ) 


C=39.4784/(A**2) 
VOL=(A**3)/4. 
C1=25. 1327/VOL 

C2A=V0A/V0L 

C2B=V0B/V0L 

C2C=V0C/V0L 

C20=V0D/V0L 

ALPHAA2=ALPHAA**2 

BETAA2=BETAA**2 

GAMMAA2=GAMf'tAA**2 

ALPHAB2=ALPHAB**2 

BETAB2=BETAB**2 

6AMMAB2=GAMMAB**2 

ALPHAC2=ALPHAC**2 

BETAC2=BETAC**2 

GAMMAC2=GAMMAC**2 

ALPHAD2=ALPHAD**2 

BETA02=BFTAD**2 

GAMMAD2=GAMMAD**2 


CALCULATE  THE  G=0  COMPONENT 


VC=0. 

VD=0. 

VA=C1*(-(ZA-QA)/(2.*BETAA2)-0A/ALPHAA2+(ZA**2)/(QA*ALPHAA2))+C2A 
VB=C1*^(ZB-QB)/(2.*6ETAB2)-QB/AI  PHA82+(Z6**2)/(QB*ALPHA62)  )+C2B 


c 

IF  (XI)  100,110,100 

100  VC=C1*(-(ZC-QC)/(2.*BF.TAC2)-QC/ALPMAC2+(ZC**2)/(QC*ALPHAC2))+C2C 
C 

no  IF  (Yl)  120,130,120 

120  VD=C1*(-(ZD-Q0)/(2.*B£TAD2}-Q0/ALPHAD2+(ZD**2)/(00*ALPHAD2))+C2D 
C 

130  VR(1 )=(1.-X1)*VA+(1.-Y1 )*VB+X1*VC+Y1*VD+VR(1 ) 

VI(1)=0. 

C 

C 

C CALCULATE  OTHER  COMPONENTS 

C 

iGsgno 

c 

DO  350  1=2,181 
C 

IGS0=IT(1 ,1)**2  + IT(2,I)**2UT(3,I)**2 
IF  (IGSO.FQ.IGSOl)  GO  TO  300 
IGS01=1GS0 
KSO=C*IGSO 
C 

VC=0. 

VD=0. 

VA=C1 * ( ZA/ ( KSQ+ ( QA* ALPHAA2 ) /ZA ) - ( ZA-OA) /KSO-OA/ ( KSQ+ ALPHAA2 ) ) 
VA=VA+C2A»£XP ( -KSO/ ( 2 . *GAMMAA2 ) ) 

C 

VB=C1»(ZB/(KSO+(OB*ALPHAB2)/ZB)-(ZB-OB)/KSO-QB/(KSO+ALPHAB2)) 

VB=VB+C2B*EXP(-KS0/(2.*GAf1MAB2)) 

c 

IF  (XI)  150,160,150 

150  VC=C1*UC/(KS0+(0C*ALPHAC2)/ZC)-(ZC-0C)/KS0-QC/(KSQ+ALPHAC2)) 
VC=VC+C2C*EXP(-KSQ/(2.*GAMMAC2)) 

C 

160  IF  (Yl)  170,300,170 

170  VD=C1*(ZD/(KSQ+(OD*ALPMAD2)/ZD)-(ZD-OD)/KSO-00/(KSO+ALPHA02)) 
VD=VD+C20*EXP(-KS0/(2.*GAMMAD2)) 

C 

C 

300  CONTINUE 
C 

VCATI0N=(1.-X1)*VA+X1*VC 

VAN10N=(1.-Y1)*VB+Y1*VD 

C 

VR(I)=(VCATION+VANION)*COSGTAU(I)^C1*CR(I)/KSO+VR(I) 

VI(I)=(VCATION-VANION)*SINGTAU(I)+C1*CI(I)/KSO+VI(I) 

C 

c 

350  CONTINUE 
C 
C 

RETURN 

END 


SUBROUTINE  ITERATE  (N.X) 

C 

C 

C THIS  SUBROUTINE  PEREORMS  N ITERATIONS  OF  THE  BAND  STRUCTURE- 

C CHARGE  DENSITY  CALCULATION. 

C X IS  A RELAXATION  FACTOR  TO  IMPROVE  CONVERGENCE. 

C ! 

c 

COMMON/CHGDEN/CR{ 1 81 ) ,C! ( 1 81 ) 

COMMON/1NFO/DUMI1Y(18),A 
COMMON/TABLE/ IT( 3, 181 ),NT( 11 .11  .11 ) 

C 

DIMENSION  CRIN(181).CIIN(181) 

C 

c 

c 

C-25. 1327*4. /A**3 

c ; 

c i 

DO  500  1TERAT=1.N  1 


DO  100  1=1.181 
CRIN(I)=CR(I) 

100  CIIN(])=CJ(J) 

C 

C 

CALL  MODELV  (CR.CI) 

C 

CALL  VALCHG 
C 

C AVERAGE  OVER  THE  LAST  ITERATION  WITH  RELAXATION  FACTOR. 

C 

Y=l.-X 

C 

DO  150  1 = 1,181 
CR(I)=X*CR(I)+Y*CRIN(I) 

150  CI(I)=X*CI{I)+Y*CIIN(I) 

C 

C GENERATE  CONSISTENCY  DATA. 

C 

DELVSQ=0. 

C 

DO  175  1=2,181 
C 

GSQ=0. 

DO  170  J=1 ,3 
170  GSq=GSQ+IT(J,I)**2 
C 

DELTAVR=(CR(I)-CRIN(n)/GSQ 

DELTAVI=(CI(I)-CIIN(I))/GSQ 

C 

DELVSQ=DELVSQ+DELTAVR**2+DELTAV1**2 

C 

175  CONTINUE 
C 
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DELV=C»SQRT(DEI  VSQ) 

PRINT  200.0ELV 

200  FORMAT  {/,*  AVERAf.E  RMS  CHANGE  IN  P0TENTIAL  = *,F6.3,»  RYDBERGS*) 


500  CONTINUE 

RETURN 

END 
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SUBROUTINE  B7SAMPL 

C THIS  SUBROUTINE  SETS  UP  A STANDARD  SET  OF  K VECTORS  AND  CALLS 

C EVAL  TO  FIND  THE  ENERGIES. 

C 

DIMENSION  AK(3) 

C0NM0N/BANDS/EK(8.40) 

COMMON /INFO/  INF0BIT(16). AMAG 1 . AMAG2 
C 

C GAMMA  POINT 

C 

AK(1)=0. 

AK(2)=0. 

AK(3)=0. 

CALL  EVAL  (11,AK,.F.) 

DO  20  1 = 1 ,8 
20  EK(I,31)=EK{I.H) 

C 

C GAMMA  TO  L 

C 

DO  30  1 = 1,10 

N=ll-1 

AK(1)-I*.05 

AK{2)=I*.05 

AK(3)=I*.05 

CALL  EVAL  (N.AK,.F.) 

30  CONTINUE 
C 

C GAMMA  TO  X 

C 

AK(1)=0. 

AK(2)=0. 

DO  40  1 = 1,10 
N=I+11 

CALL  EVAL  (N,AK,.F.) 

40  CONTINUE 
C 

C X TO  K 

C 

AK(1)=1, 

AK(2)=.125 

AK(3)=.125 

CALL  EVAL  (22,AK,.F.) 

C 

C < TO  GAMMA 

C 

AK(3)=0. 

DO  50  1=1,8 
N=31-I 

AK(1  ) = . 09375*1 
AK{2)=, 09375*1 
CALL  EVAL  (N,AK,.F.) 

50  CONTINUE 
C 
C 
C 
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RANGE =20. 

^ AZER0=5.*AINT(EK(3,11)/5.)-10, 

^ CALL  BANDPLT  (AZERO, RANGE) 

C 

RETURN 


C 

END 
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SUBHOUTINE  PR  INGAP 

■ C 

; C 

’ C THIS  SlIRROUTINE  CALCULATES  THE  PRINCIPAL  BAND  GAPS  AT  GAMMA, 

f C L AND  X. 

C 

[ C 

t DIMENSION  AK(3) 

[ DIMENSION  GAP(8) 

C0MV0N/8ANDS/E(8,40) 

' COMMON/HEADING/NAME (8). I DATA (8) 

REAL  L3V,L1C,L3C 
C 
C 

DO  20  I=I .3 
20  AK(I)=0. 

C 

CALL  EVAL  (n.AK,.E.) 

C 

AK(3)=I. 

C 

CALL  EVAL  (21,AK,.F.) 

C 

DO  30  1=1,3 
30  AK(I)=.S 
C 

CALL  EVAL  (1,AK,.F.) 

C 
C 

ENTRY  GAPS 
C 
C 

C SORT  OUT  CiAMMA  LEVELS 

C 

G15V=E(3,n) 

G1C=E(5,11) 

G15C=E(6,n) 

C 

IF(ABS(E(2,n)-G15V).LE.0.01)  GO  TO  100 
C 

G1C=E(2,11) 

PRINT  80 

80  FORMAT  {*  METAL  BAND  STRUCTURE*) 

C 

100  IF  (ABS(G1C-G15C).LE.0.01)  G1C=E(8,11) 

C 

C L and  X LEVELS 

C 

L3V=E(4,1) 

LlC=E(b,l 
L3C=E(6,1) 

C 


C 


X4V=F.(4,21) 

X1C=E(5,21) 

X3C=E(6,21) 


CALCULATE  GAPS 


GAP(1)=G1C-G15V 

GAP(2)=G1SC-G15V 

GAP(3)=L1C-G15V 

GAP(4)=X1C-G]5V 

GAP(5)=L1C-L3V 

GAP(6)=L3C-L3V 

GAP{7)-X1C-X4V 

GAP(8)=X3C-X1C 


PRINT  200 

PRINT  210,(GAP(I).I=1 ,fi) 

PRINT  220,{inATA(I),I=I,8) 

200  FORMAT  (//,12X,*G1C-G15V  G15C-G15V  L1C-G15V  XIC-G15V  L1C-L3V 
1 L3C-L3V  X1C-X4V  X3C-X1C*,/) 


210  FORMAT  (*  RESULTS  *,8F10.2) 
220  FORMAT  (/,»  DATA  *,8A10) 


oooo  oo  o o oo  o o o o o o o o oo  o ooooooo 


SUBROUTINE  VALCHG 


THIS  SUBROUTINE  EVALUATES  THE  VALENCE  CHARGE  DISTRIBUTION  USING 
THE  CHADI-COHEN  3 POINT  BRILLOUIN  ZONE  SAMPLE. 


COMMON/CHGDEN/CR( 1 81 ) ,CI ( 1 81 ) 
C0MM0N/INF0/INF0BIT(16) 
C0MM0N/TA2LE/IT(3,181).NT(11  ,11  ,11 ) 
LOGICAL  INFOBIT 

DIMENSION  AK(3) 


DO  100  1 = 1,181 
CR(1)=0. 

100  CI(I)=0. 

AK(1)=.5 

AK(2)=0. 

AK(3)=0, 

CALL  EVAL  (33,AK,.T.) 
CALL  ABSPS12  (.25) 

AK(2)=.5 

CALL  EVAL{34,AK,.T.) 
CALL  ABSPSI2  (.5) 

AK(1)=1. 

CALL  EVAL  (35,AK,.T.) 
CALL  ABSPSI2  (.25) 

CALL  SYMETRZ 


PRINT  CHARGE  COEFFICIENTS  IF  INF0BIT(4)  IS  TRUE 

IF  (.N0T.INF08IT(4))  GO  TO  300 
C 

PRINT  200 
PRINT  205 

PRINT  210,(I,(IT(J,I),J=1,3),CR(I),CI(I),I=1,181) 


200  rORMAT  (///,20X,*F0URIER  COMPONENTS  OF  THE  VALENCE  CHARGE  DENSITY 


o o o 


107 


1.//) 

205  FORMAT  (4(*  MO  ( I J K)‘ , 7X  .»CR* ,6X  ,*CI* ,4X) ./) 
210  FORMAT  (4(14/  (*C3I2,*)  *.2F8.4,4X)) 

C 

300  CONTINUE 

PLOT  CHARGE  DENSITY  IF  1NF0B1T(5)  IS  TRUE 

IF  (1NF061T(5))  CALL  FSPLT(0. ,40. ) 

C 

C 

RETURN 

END 
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SUBROUTINE  EVAL  (NO. AK. CONTROL) 

THIS  SUBROUTINE  EVALUATES  THE  BAND  ENEPOIFS  FOR  A GIVEN  K USING 
THE  POTENTIAL  PART  OF  THE  HAM;'.  TOMAN  CALCULATED  PREVIOUSLY  AND 
STORED  IN  POTNTL. 

N IS  THE  INDEX  UNDER  WHICH  THE  RESULTS  ARE  STORED  IN  'BANDS' 

AK  IS  THE  K VECTOR  AT  WHICH  THE  ENERGIES  ARE  EVALUATED. 

A IS  THE  LATTICE  CONSTANT. 


LOGICAL  CONTROL 
LOGICAL  INFOBIT 
LOGICAL  INC 

DIMENSION  AK(3) 

DIMENSION  T(181),AKSQ(181),INDX(18i) 

DIMENSION  D(A0),E l(40),E2(A0).TAU(2,40) 

DIMENSION  ZR(1600),2I(I600} 

DIMENSION  AR(181 .8).AI(181 ,8),BR(181 ,8), 81(181 ,8) 

COFWON/ INFO/ INFOB  I T( 16 ) .AMAGl .AMAG2 ,A 
C0MM0N/BANDS/E.<(8,40) 

COMV'lON/LOCALV/VRdai ) ,VI  (181 ) 

COMMON/TABLE/ ITAB(3 ,181) ,NTAS (11,11,11) 

COMMON/HAMLT/HR( 1600) ,HI ( 1600) 

COMMON//VMNR(  1500)  .VMM  ( 1600) , VJMR(3200) , VJMl  (3200) , VJJ(80) 
C0MM0N/INCLUD£/INC(181) 

EQUIVALENCE  (VMNR.AR) ,(VMNI ,AI) , (HR,BR) .(HI ,BI ) 

EQUIVALENCE  ( VHN R , Z R ) , ( VMN I . Z I ) 

EQUIVALENCE  (T.AKSQ) 


DO  20  1 = 1 ,181 
AKSQ(I)=0. 

INDX(I)=I 
00  20  J=l,3 

20  AKSQ(I)=AKSQ(I)+(AK(J)+ITAB(J,I))**2 


CALL  SORT  (181 .AKSQ.INDX) 

C0NST=(6.2832/A)**2 
DO  25  1 = 1 ,181 
IF(AKSQ(I).LE. AMAGl)  N1=I 
1F(AKSQ(1).LE.AMAG2)  MAX=I 
25  T(I)=C0NST*AKSQ(I) 

N2=MAX-NI 

IF  (N1.LE.40)  GO  TO  35 
PRINT  30. NO 

30  FORMAT  (*  INNERSET  OVERRANGE  ON  *,I5) 


oooo  o o o ooo  ooo  ooooooo 


Nl=40 

N2=MAX-N1 

c 

35  IF  (N2.LE. 3200/Ml)  GO  TO  40 
C 

PRINT  37. MO 

37  F0P.MAT  (*  OUTER  SET  OVERRANGE  ON  *.I5) 
N2=3200/N1 
C 

40  CONTINUE 
C 
C 

IF(INF0BIT(1))  CALL  IMTYME 


SET  UP  HAMILTONIAN  FOR  L0WDIN3 


SMALL  SQUARE  MATRIX 

DO  60  N=1,N1 
DO  60  M=1  ,N 
NM=N1*(M-1  )+N 
MN=N1*(N-1 )+M 

CALL  POTNTL ( 1 NUX { M) , I NDX ( N ) , VMNR ( MN ) , VMN I (MN ) ) 

VM,NR(NM)=VMNR(MN) 

VMNI(NM)=-VMNI(MN) 

60  CONTINUE 

ADD  KINETIC  ENERGY  TO  DIAGONAL  ELEMENTS 

DO  65  N=1  ,N1 
NN=N1*(N-1)+N 
65  VMNR(NN)=VMNR(NN)+T(N) 

DIAGONAL  ELEMENTS  OF  LARGE  SQUARE  MATRIX 

DO  70  J=1 ,N2 
70  VJJ(J)=VR(1)+T(N1+J) 

RECTANGULAR  MATRIX 

DO  80  J=1 ,N2 
DO  80  M-1  ,N1 
n=iNnx(j+Ni) 

I2=INDX(M) 

JM=N2*(M-1)+J 

CALL  POTNTL (II, I2,VJMR( JM) , VJMI ( JM) ) 

80  CONTINUE 

EBAR  IS  THE  AVERAGE  OF  THE  FIRST  8 ENERGY  LEVELS  OF  THE  EMPTY 
LATTICE. 

EBAR=0. 

DO  85  J=1 ,8 

85  EBAR=EBAR+VMNR(N1*J-Nl+J)/8. 


no 


c 

ir(INFOBIT(D)  CALL  TYf'ER(5HSETUP) 

CALL  LOUniN3(rn  .NC.FBAfi) 
lF(INFOniT(l ))  CALL  TYMER{7HLOWDirj3) 

C 

C 

C CALL  DIAGONALIZATION  ROUTINES 

C 

IF  (CONTROL)  GO  TO  300 
C 
C 

CALL  HTRIDI  (Ml  ,M1  .HR,[|I  .D.El  ,E?,TAU) 

IF  (INFOBIT(D)  CALL  TYMER  (6HHTRIDI) 

C 

CALL  IMTQLl (N1 .D.El ,IERR) 

IF  (INFOBIT(I))  CALL  TYIIER  (6HIMTQL1 ) 

C 

C 

C IF  ERROR  CODE  INDICATES  EIGENVALUE  DID  NOT  CONVERGE,  PRINT  OUT 

C MESSAGE  AND  TERMINATE  EXECUTION. 

C 

IF  (lERR.NE.O)  GO  TO  900 
C 
C 

C CONVERT  ENERGIES  TO  EV  AND  STORE  IN  OUTPUT  ARRAY 

C 

140  DO  150  1=1 .8 

150  EK{I,N0)=13.605*D(I) 

C 

IF  (.N0T.INF0BIT(3))  GO  TO  190 
C 

PRINT  170.N0,(AK(J),J=1 ,3) 

PRINT  171,N1,N2 

PRINT  172,(EK(I ,N0),I=1 ,8) 

170  FORMAT(//,*  N0.*,I5,10X,*K=  *.3F5.2) 

171  F0RMAT(I5,*  PLANE  WAVES  TREATED  DIRECTLY,*, 15,*  PLANE  WAVES  TREATE 

ID  BY  PERTURBATION*) 

172  FORMAT  (/,*  ENERGIES  *,8F10.2) 

190  CONTINUE 

C 

RETURN 

C 

C 

C EIGENVECTOR  CALCULATION 

C 

C INITIALIZE  Z 

Q 

300  DO  305  1=1,1600 
ZR(I)=0. 

305  ZI(I)=0. 

C 

DO  310  1 = 1 ,N1 
310  ZR(I+N1*I-N1)=1. 


C 


IF  (INFOBIT(l))  CALL  INTYME 


1 i i ) 


c 

C CALL  DIAGONAL IZATION  ROUTINES 

C 

CALL  HTRIDI  (M  ,N1  .HR.HI.D.El  .F2.TAU) 
IF  (INFOBIT(l))  CALL  TYMER  (GitHTRIDI) 
C 

CALL  IMTQL2  (N1  .NT ,D,E1 .ZR.IERR) 

IF  (INFOBIT(I))  CALL  TYMER  (6HIMTQL2) 
C 

CALL  HTRIBK  (N1  ,N1 .HR.HI ,TAU,4 ,ZR,ZI ) 
IE  (INFOBIT(l))  CALL  TYMER  (6HHTRIBK) 
C 
C 

C TRANSFORM  AND  OUTPUT  EIGENVECTORS 

C 

C 

NBAND--4 

IF  (INFOBIT(16))  NBAND=8 
NWAVE=181*NBAND 
C 

DO  315  I = 1,N»1AVE 
C 

BR(n=0. 

316  BI(I)-0. 

C 

C 

MDISP=0 

NOISP=0 

C 

C 

DO  325  IBAND=1,NBAND 
C 

DO  320  1 = 1 , N1 
C 
C 

BR(I+MDISP)=ZK(I+NDISP) 

320  3I(HHDISP)=ZI{I+NDISP) 

C 

MDISP=HDISP+181 
325  NDISP=NDISP+N1 
C 
C 

C UNFOLD  THE  LOWDIN  TRANSFORMATION. 

C 

DO  350  J=1,N2 
C 

DEN0M=1./(EBAR-VJJ(J)) 

MDISP=0 

C 

C 

DO  345  IBAND=1  ,NBAND 
C 

CR=0. 

CI=0. 

JM=J 

C 


oo  o o o 


DO  340  KDEX=1 .Nl 
C 

M=M01SP+MDEX 

C 

CR=CR+VJMR{  JM)*BR(M)-VJMI  (f1) 

CI=CI+VJMI{JM)*BR(«)+VJnR(JM)*BI(M) 
C 

340  JM=JM+N2 
C 
C 

BR(j+Ni+rinisn^-CR*nENnn 

BI(J»N1+MDISR)=C:*DEN0M 

c 

345  MDISP=riDISP+181 
C 

350  CONTINUE 


REARRANGE  COMPONENTS. 


DO  380  1=1,181 
C 

J=INDX(I) 

C 

INC(J)=I,GT,MAX 

ID=I 

C 

C 

DO  380  IBAND=1  .NBAND 
C 

AR(J)=BR(1D) 

AI(J)=BI(ID) 

C 

J=J+181 

ID=ID+181 

C 

380  CONTINUE 
C 
C 

GO  TO  140 
C 

900  PRINT  910.no 

910  FORMAT  (//,*  DIAGONAL IZATION  ROUTINE  FAILED  ON  *.15,*  K VECTOR*) 

999  FORMAT  (/////) 

END 


oooo  ooooooo 


SUBROUTINE  P0TMTL(M,N.V1R,V1 1 ) 


li'i 


THIS  SUBROUTINE  CALCULATES  ”HE  M.N  ELEMENT  OE  THE  POTENTIAL 
MATRIX  ANO  STORES  THE  RESUIT  IN  VI.  NOTE  THAT  M,N  HERE  REFER 
TO  INDICIES  AS  ORIGINALLY  DEFINED  IN  'TABLE'. 


COMMON/', -OCALV/VR(  181 ) ,VI  (181 ) 
COMMON/TABLE/ I (3, 1 81 ),NTAB(11 ,11  ,11) 


INDEX  MANIPULATIO.NS 


L1=I(1,M)-I(1,N) 

L2=l(2,M)-I(2,N) 

L3=I(3,M)-1(3,N) 

IF  ((L1**2+L2**2+L3**2).GT.2A)  GO  TO  20 
C 

10  NA=NTAB(L1+G,L2+6,L3*6) 

V1R=VR(NA) 

V11=VI(NA) 

GO  TO  30 
C 
C 

20  V1R=0. 

V1I=0, 

30  CONTINUE 
C 

RETURN 

END 


^ OOO  O OOOOOO  O O o OOOOOOO 
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SUBROUTINE  L0W0IN3  (N1,N2,EBAR) 


N1  IS  THE  DIMENSION  OF  THE  SMALL  SQUARE  MATRIX  AND  N2  IS  THE 
DIMENSION  OF  THE  RECTANGULAR  MATRIX.  N1  MUST  BE  LE  40  AND 
N2  mST  BE  LE  80. 


DIMENSION  DENOM  (80) 

COMMON// VMN  R { 1 600 ) . VMN I ( 1 600 ) , V JMR ( 3200 ) , V JMI ( 32D0 ) , VJ J ( 80 ) 
C0MM0N/t6U1LT/HR(1600).HI(1600) 


00  25  0=1, N2 
DENOM(J)=EBAR-VJJ(J) 
25  CONTINUE 


CALCULATE  THE  OFF  DIAGONAL  PERTURBATION  TERMS 

TRIANGULAR  NESTED  DO  LOOPS 

DO  40  M=2,N1 
MAX=M-1 
INDXM=N2*MAX 
DO  40  N=1,MAX 
INDXN=N2*(N-1 ) 

PERTR=0. 

PERTI=0. 

DO  30  0=1 ,N2 

M0=0+INDXM 

N0=0+INDXN 

PERTR=PERTR+(VOMR(MO)*VOMR(NO)+VOMI(MO)*VOMI(NO))/DENOM(0) 
PERTI=PERTI+(VOMR(MO)*VOMI(NO)-VOMI(MO)*VOMR(NO))/DENOM(0) 
30  CONTINUE 

ADO  PERTURBATION  TERM  TO  HAMILTONIAN 

MN=N1*(N-1)+M 

NM=N1*(M-1)+N 

HR(MN)=VMNR(MN)+PERTR 

HR(NM)=HR(MN) 

H1(MN)=VMM(MM)+PERTI 

HI(NM)=-HI(MN) 

40  CONTINUE 

CALCULATE  THE  DIAGONAL  MATRIX  ELEMENTS 

DO  60  N=1 ,N1 
E=VMNR(N1*(N-1)+N) 

PERT=0. 

DO  50  0=1  ,N2 
0N=N2*(N-l)+0 

PERT=PERT+(VOMR(ON)**2+VOMI(ON)-‘*2)/(E-VOO(0)) 


50  CONTINUE 
NN=NI»N-NWf( 
HR(NN)=VMNR(NN)+PERT 
HI(NN)=0. 

60  CONTINUE 
C 
C 

RETURN 

END 


,1 


1 
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SUBROUiINC  liANDPLT  {AZCRO.RANGl  ) 

C 

C THIS  SUBROUTINi;  GENCRATES  A PRINTER  PLOT  OF  THE  BAND  STRUCTURE 

C 

C 

C0W»N/BANDS/EK(8,40) 

COFWON/HEAOING/NAML(8) 

COMMON// ICHAR( 31  ,51 ) ,ESCALE(51 ) 

C 

DELTA=P,ANGE/50. 

DO  25  1=1,1581 
25  ICHAR(I)=0. 

C 

C SET  UP  EACH  COLUMN 

C 

DO  50  N=l,31 
C 

C SCAN  UP  THE  COLUMN 

C 

00  50  11.51 

AMIN  •VER0»(I-1.5)*0ELTA 
AMAX  AZERnt(I-.5)*DELiA 
C 

DO  50  J = l,9 
X=LK(J,N) 

IF((X.LE.AMAX).AND.(X.GT.AMIN))  ICHAR( N , I ) = ICHAR( N , I )+l 
50  CONTINUE 
C 

DO  55  1 = 1,51 

55  ESCAL£(1)=AZER0+(I-1)*DELTA 
C 

C CONVERT  INTEGERS  TO  CHARACTERS 

C 

DO  70  1 = 1 ,1581 
ICHAR(I)=ICHAR(l)+0000033 
IF( ICHAR{ I ) .EQ. 0000033)  ICHAR( I )=0000055 
70  CONTINUE 
C 

C GENERATE  OUTPUT 

C 

PRINT  100,(NAME(I),I=1 ,8) 

100  F0RI-1AT  (*1BAND  STRUCTURE  FOR  *,8A10,//) 

1=51 

DO  150  NA=1,51 
N=5?-NA 

IF(N.EO.I)  GO  TO  120 
PRINT  200,(ICHAR(J,N),J=1,31) 

GO  TO  150 

120  PRINT  201,ESCALE(I),(ICHAR(J,N),J  = 1 ,31) 

1=1-5 

150  CONTINUE 
C 

200  F0RMAT(13X,*1*.R1 ,30(3X,R1 )) 

201  F0RMAT(2X,n0.2,*  I*,R1 ,30(3X,R1 )) 


C 


PRINT  210 
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PRINT  212 
PRINT  211 


210  F0R/WT(1<}X.*I*,39(*,*).*I*.39(*.»).*I*,7(*.*),*I*,3U*  *)  *I*) 

211  F0RMAT(UX.*L*,17X,*(111  )*,15X.*GAHMA*,15X,*(100)*,17X  *X*,7X 
1*K*.11X, *(110)*, 12X, ‘GAMMA*) 

212  FORMAT (/) 


RETURN 

END 


SUBROUTINE  AGSPSI2  (ALFA) 


THIS  SUBROUTINE  CALCULATES  THE  ABSOLUTE  SQUARE  OF  A BLOCH 
FUNCTION. 


LOGICAL  INC 

C0MM0N/TABLE/IT(3,181).NT(11 .11  ,11) 

COMMON  /HAMLT/HR(1600),HI{lo00) 

COMMON/ /VtlNP, ( 1 600 ),  VMN I ( 1 600 ) 

COMMON/CHGDEN/CR( 181 ) ,CI  ( 181 ) 
COMMON/INCLUDE/INCdCl) 

DIMENSION  AR(181,8).AI{181.8),BR(181 ,8) ,BI ( 181 ,8) 
EQUIVALENCE  (VMNR.AR) , (VMNI ,AI ) , (HR.BR) . (HI ,BI ) 


DO  80  1=1  ,724 

BR(I)=0. 

BI(1)=0, 


ID=0 

DO  210  1=1,181 
IF  (INC(D)  GO  TO  210 


JD=0 

DO  200  J=l,181 
IF  (INC(J))  GO  TO  200 


I1=IT(1+ID)-IT(1+J0) 

I2=1T(2+ID)-1T(2+JD) 

I3=IT(3+ID)-IT(3+JD) 

I5Q=I1**2+I2**2+I3**2 

IF  (ISQ-32)  100,100,200 

M=NT(Il+6, 12+6, 13+6) 


DO  150  IBAND=1,4 


BR(M)=BR(M)+AR(11 )*AR(J1 )+AI (1 1 )*AI (J1 ) 
BI{M)=B1(M)+AI(I1)+AR(J1)-AR{I1)*AI(J1 ) 

M=M+181 
11  = 11+181 


I 


lis 


150 


200 

210 


J1=J1+181 

CONTINUE 

TD=JD+3 

10=10+3 


NORMALIZE  TO  1 ELECTRON  PER  STATE. 
1 0=0 

DO  250  IBAND=1,4 

C0NST=ALEA/BR(1 ,IB4ND) 

DO  240  1 = 1 ,181 

CR(I)=CR(I)+CONST*BR(I+ID) 
CI(I)=CI(I)+CONST*BI(I+ID) 

240  CONTINUE 

250  10=10+181 


I 


RETURN 

END 


o O Oo  oooo  o o oooo  o o oooooo 


SUBROUTINE  SYMETR2 


THIS  SUBROUTINE  AVERAGES  THE  CHARGE  DENSITY  OVER  THE  POINT 
GROUP  OPERATIONS. 


C0MM0N/TABLE/IT(3,181),NT(11 ,11 ,11  ) 
COMMON/SFACTR/COSGTAU{ 181 ) ,SINGTAU( 181 ) 
COMMON/CHGDEN/CR( 181 ) ,CI (1 81 ) 

COMMON/ PTGROUP/  Gf1N(3,3,24) 

COMMON //AR ( 200 ) , A I ( 200 ) , BR { 200 ) , B I ( 200 ) 

DIMENSION  L(3) 

EQUIVALENCE  (L( 1 ) ,L 1 ) , (L(2)  ,L2) , (L (3)  ,L3) 


TRANSLATE  THE  COORDINATE  SYSTE*'  TO  A SYMMETRY  POINT. 
DO  100  1 = 1,181 

BR(  1 ) = CR( I )*C0SGTAU{ I ) -Cl ( I )*SINGTAU( I ) 
BI(I)=CR(I)*SINGTAU(I)+C1(I)*C0SGTAU(I) 

AR(I)=0. 

100  AI(I)=0. 


AVERAGE  OVER  THE  GROUP  OPERATIONS. 


JD= 

=0 

DO 

250 

0=1 

.24 

ID= 

=0 

DO 

200 

1 = 1 

,181 

KD= 

=0 

00 

180 

K=1 

.3 

L(K)=6 

DO 

150 

M=1 

.3 

150  L{K)=L(K)+GMN(M+KD+JD)*IT(M+ID) 
C 

180  KD=KD+3 
C 
C 

N=NT(L1,L2,L3) 

C 

AR(I)=AR(I)+BR(N)/12. 

AI(I)=AI(I)+BI(N)/12. 


o o o o 
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c 

200  1D=I0>3 
C 

250  J0=JD+9 


TRANSLATE  THE  COORDINATE  SYSTEM  BACK. 

DO  300  1=1  ,181 
C 

C R ( I ) = AR  ( I ) *COS  r.TAU  ( I ) + A I ( I ) *S  I NOT  AU  ( I ) 
300  Cl  { I ) =AI  ( I ) *C0Sr,TAU(  I ) -AR(  I )*SI  NGTAU(  I ) 

C 

C 

RETURN 

END 


SUBROUTINE  SLATER3  (ALPHA, CR, Cl .XSR.XSI ] 


THIS  SUBROUTINE  CALCULATES  A QUADRATIC  APPROXIMATION  TO  THE 
SLATER  EXCHANGE  APPROXI'IATiON. 

THE  INPUT  MUST  BE  IN  UNITS  OF  ELECTRONS  PER  PRIMITIVE  UNIT  CELL 
AND  THE  OUTPUT  IS  IN  RYDBERGS. 


DIMENSION  XSR(18)).XSin81),CR(18'),CI(181) 

C0MM0N/INF0/DUMMY(18),A 

C0=0.8A7 

C1=0.M2 

C2=-0. 00209 

C0NST=-4.68956*ALPHA/A 


CALL  CONVOL  (181  ,CR,CI .CR.CI .XSR.XSI ] 


DO  100  1 = 1,181 

XSR(  I ) = C0NST*(C1*-CR(  I )+C2*XSR{  I ) ) 
XSI(1)-C0NST*(C1*CI(I)+C2*XSI(I)) 

CONTINUE 

XSR(1)=XSR(1)^C0NST*C0 


RETURN 


SUBROUTINE  CONVOL  (N.AR.AI .BR.BI  .CR.r  : ) 


C 

C 

C 

C THIS  SUBROUTINE  CALCULATES  THE  FOURIER  SERIES  FOR  THE  PRODUCT 

C OF  TWO  FOURIER  SERIES. 

C 

C AR  AND  AI  ARE  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  FIPST  SERIES, 

C BR  AND  BI  ARE  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  SECOND  SERIES 

C CR  AND  Cl  ARE  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  PRODUCT 

C SERIES. 

C 

C 

DIMENSION  AR(M),AI(N),BR(U),BI{N),CR(N).CI(N) 

C 

COMMON/TABLE/ITO.ISI  ),NT(11  ,11,11) 

C 

C 

MAX=0 

DO  no  1 = 1,3 
no  MAX=MAX+IT(I  ,N)**2 

c 

c 

DO  1?0  1 = 1, N 
CR(I)=0. 

120  CI(1)=0. 

C 

C 

DO  200  1=1 ,N 
DO  200  J=1 ,N 
C 
C 

ii=iT(i,i)nT(i,j) 

I2=1T(2,I)+IT(2,J) 

I3=IT(3,I)+IT(3,J) 

C 

130=11*11^12*12+13*13 
IF  (ISQ.GT.MAX)  GO  TO  200 
C 

M=NT(Il+6, 12+6, 13+6) 

C 

CR(M)=CR(M)-iAR(I)*BR(J)-AI(l)*BI- J) 

C1(M)=CI(M)+AR(I)*BI (J)+AI(1)*BR(J) 

C 

200  CONTINUE 
C 

c 

RETURN 

END 
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SUBROUTINE  FSPLT  (ZERO .RANGE ) 

C 

C 

C CR  IS  THE  ARRAY  OF  THE  REAl  COEFICIENTS. 

C Cl  IS  THE  ARRAY  OF  THE  If'AGINARY  COEFICIENTS. 

C A IS  THE  STARTING  POINT  IN  THE  UNIT  CELL. 

C B IS  THE  FINAL  POINT. 

C 

C 

C0KM0N/CH6DEN/CRn81),CI(131 ) 
C0Mf'10N/TABLE/I{3,181  ),NTAB(11  ,11 ,11 ) 
C0MM0N/HEADiNG/NAME(8) 
C0MM0N//X(101),Y(101),ICHAR(101) 

C 

DIMENSION  A(3),B(3),DELX(3) 

C 

C 

N = 181 
C 

DO  no  J=1  ,3 
A(J)=0. 

B(J)=1. 

110  0ELX(J)-(B(J)-A(J))/100, 

C 

DO  115  L=l,101 
115  Y(L)=0. 

C 

C 

C 

DO  200  K=1,N 
G0X=0. 

C 

DO  130  J=l,3 

130  GDX=GDX+I(J,K)*DELX(J) 

C 

CSGDX=C0S(6.2832*GDX) 

SNGDX=SIN(6.2832*G0X) 

C 

C 

C CALCULATE  TERM  AT  A 

C 

GX=0. 

DO  135  J=l,3 
135  GX=GX+I(J,K)*A(J) 

C 

CSGX=C0S(6.2S32*GX) 

SNGX=SIN(6.2832*GX) 

C 

Y{1)=Y(1)+CR(K)*CSGX-C!(K)*SNGa 

C 

C CALCULATE  TERMS  FROM  A TO  B 

C 

DO  200  L=2,ini 
C 

CSTEMP=CSGX 

SNTEMP=SNGX 
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CSGX=CSTEMP*CSGDX-SNTEMP*S‘.'GDX 

SNGX=SNTEMP*CSGOX+CSiEMP*SNGOX 

C 

Y(L)=Y(L)+CR(K)*CSGX-CI{K)*StU  X 
C 
C 

200  CONTINUE 
C 
C 
C 

C GENERATE  OUTPUT 

C 

PRINT  205,(NAriE(J),J  = l,8) 

205  FORMAT  (*IPLOT  OF  *,8A10,/) 

PRINT  206,(A(J),J=1 ,3),(B(J),J=1 ,3) 

206  FORMAT(*  FROM  (*,F6.3,*,*,F6.3,*.*,F6.3,»)  TO  (*,F6.3,*,*, 
1F6.3,*,*,F6.3,*)*,//) 

C 

C 

0ELY=RANGE/50. 

1FLG=50 

C 

DO  300  1N0X=1,51 
C 

M=51-INDX 

YM=M*DELY+ZERO 

YMIN=YM-.5*DELY 

YMAX=YM+.5*DELY 

C 

C 

DO  220  K=l,101 
C 

ICHAR(K)=0000055 

IF  ((Y{K).LE,YMAX).AND.(Y(K).GT,YMIN))  ICHAR(k)=0000047 
C 

220  CONTINUE 
C 

c 

c 

IF  (M.EQ.IFLG)  GO  TO  230 
C 

PRINT  225,(ICHAR(K),K=1  ,101) 

225  FORMAT  (15X,*I*,101R1 ) 

GO  TO  300 
C 
C 

230  PRINT  235.Yn.(ICHAR(K),K=l ,101) 

235  FORMAT  (5X,F10.2,*I*,101R1 ) 

IFLG=lFLG-5 

C 

300  CONTINUE 
C 
C 

PRINT  310 

310  FORMAT  (15X ,20(*I ....*) .*!*) 
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PRINT  311  ,(J,J  = 1 ,9) 

311  FORMAT  (13X,‘0.0*,9(7X,*0.*.]1).7X.*1.0*) 
C 

c 

DO  350  L=1 ,101 
350  X(L)={L-1)/100. 

C 

PRINT  999 

PRINT  355,(X(L),Y(L),L=1,101) 

355  FORMAT  ( 1 X ,5{F5. 2,F10. 3, lOX) ) 

C 

999  FORMAT  {/////) 

C 

C 

RETURN 

END 
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SUBROUTINE  TYMER  (NAME) 

C 

C THIS  SUBROUTINE  GENERATES  TIMING  INFORMATION 

C 

C 

CALL  CPELAP(I) 

T={lMILLI(I)-K)/IOOO. 

PRINT  lO.T.NAME 

10  FORMAT  (5X,F10.3,*  SECONDS  REQUIRED  TO  EXECUTE  *,A10) 
C 

ENTRY  INTYME 
CALL  CPELAP(I) 

K=IM1LLI(I) 

RETURN 

END 


SUBROUTINE  SORT(N,X,I) 

C 

C THIS  SUBROUTINE  SORTS  A ONE-DIMENSIONAL  ARRAY  X IN  ORDER  OF 

C INCREASING  MAGNITUDE  AND  SIMULTANEOUSLY  PERMUTES  AN  INDEX  ARRAY 

C I.  N IS  THE  LENGTH  OF  THE  ARRAY 

C 
C 

DIMENSION  X(N),I(N) 

DO  50  J=2,N 
XT=X(J) 

IT=I(J) 

DO  40  K=2,J 
L=J-K+1 

IF(XT,GE.X(D)  GO  TO  50 
X(L+1)=X(L) 

I(L+1)=I(L) 

X(L)=XT 
I(L)=IT 
40  CONTINUE 
50  CONTINUE 
RETURN 
END 
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PROGIt/'.M  CRVSPOT  ( INPUT  .OUTPUT  .TAPCl  = INPUT  ,TAP[6) 

C 

c 

c 

COMMON/TABLE/ IT(3,Un  ),NT(  11  .n  ,11 ) 

COMMON/PARAf.TR/ZA.gA.ALPHAA.BETAA.GAI'.MAA.VOA.ZB.QB.ALPHAB.BETAB, 
IGAMmB.VOB.XAUPHA.HELAX.El  ,E2,A 
C0MM0N/SFACTR/C0SGTAU(181  ),Sir.GTAU{181 ) 

COMMON/DI FF/AR( 181  ),A 1(181) 

COr-WON/CHGDCN/CR  { 1 81 ) ,C  I ( 1 81 ) 

C0Mn0N/HEADING/f.'AME(8) , IDATA(8) 

C0MM0N/INF0/INF(16) 

COMMON/PRI1TR2/ZC,QC,ALPHAC,BETAC,GAriMAC,VOC,ZD.QD,ALPHAD,EETAD, 
IGAKWD.VOD.Xl  ,Y1 
C 

DIMENSION  XA(3),Xn(3).X(3) 

DIMENSION  PP,MTR(17l 
DIMENSION  PRHT2(14) 

C 

EQUIVALENCE  (PRMTR( 1 ) ,ZA) 

EQUIVALENCE  (PRnT2( 1 ) ,ZC) 

C 

C 

logical  inf 
DO  10  1=1,16 
10  INF=.F. 

C 

C 

C 

READ  60.(NAME(I),I=1,8) 

IF  (E0F,1)  100,50 
50  READ  60,(IDATA(I) ,1=1 ,8) 

READ  80,(PRMTR(I),I  = 1 ,17) 

READ  90,(J.CR(J).CI(J),I=1,181) 

C 

60  FORMAT  (8A10) 

80  FORMAT  (8F10.4) 

90  FORMAT  (4(14, 2F8. 5)) 

C 

DO  95  1 = 1 ,14 
95  PRMT2(1)=0. 

C 

GO  TO  150 
C 

100  CALL  GET  (5KTAPE6,7HPARAMTR,0,0) 

C 

READ  (6,110)  (PRMTR(I),I=1 ,17) 

READ  6,120)  (NAME(I),I=1,8),(1DATA(I),I=1,8) 

READ  (6,110)  (PRMT2(IKI  = 1,14) 

I REWIND  6 

’ C 

CALL  GET  (5HTAPE6,7HCHGDATA,0,0) 

READ  (6,130)  (CR(I),Cl(I),I=i,181) 

C 

I 

! 
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no  FORMAT  (FI 2. A) 

120  FORMAT  (8A10) 

130  FORMAT  (2F10.5) 

C 

c 

150  CONTINUE 
C 

CALL  FCCLAT 
C 
C 

CALL  DATE  (11) 

CALL  TYME(I2) 

PRINT  170,11,12 

PRINT  175,(NAME(I),I=1 ,8) 

PRINT  180,ZA,QA 

PRINT  185, Al PHAA,BETAA,GAMMAA,V0A 
PRINT  190,7B,QB 

PRINT  185,  ALPHAS,  BETAB,GAM:IAB,V0B 
PRINT  195,X1,Y1 
PRINT  999 
C 

170  FORMAT  (10X,*RUN  DATE  *,A10,10X,*TIME  *,A10,//) 

175  FORMAT  (*  ELECTROSTATIC  POTENTIAL  OF  *,0A10,//) 

180  FORMAT  (*  A ATOM*,/,*  Z=*,F2.0,5X,*Q=*,F2.0) 

185  FORMAT  (*  ALPHA=*,F5.2,5X,*BETA=*,F5.2,5X,*GAMMA=*,F5.2,5X,*V0=*,F6. 

16.2) 

190  FOR.MAT  (//,*  B ATOM*,/,*  Z-*,F2.0,5X ,*Q=*,F2.0) 

195  FORMAT  (/,5X,*X=*,F4.2,10X,*Y=*,F4.2) 

C 

C 

C 

CALL  FOURCF 
C 

DO  350  J=1 ,3 
350  X(J)=.375 

CALL  CLUSTER  (A.X , 10. ,X1  ,Y1 ,VREF1 ) 

CALL  FSPT  (181 ,X,AR,AI ,VREF1 ) 

C 

00  360  0=1,3 
360  X(J)=.625 

CALL  CLUSTER  (A,X,10.  ,X1  ,Y1 ,VREF2) 

CALL  FSPT  (181  ,X,AR,AI,VREF2) 

C 

VREF1=-13.605*VREF1 

VREF2=-13.605*VREF2 

C 

PRINT  375.VREF1  .VREF2 

375  FORMAT  (*  VREF1=*,F10.2,10X,*VREF2=*,F10.2) 

C 

C 

900  CONTINUE 
999  FORMA!  (/////) 

C 

END 


} 
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SUBROUTINE  FOURCF 
C 

c 

c 

common/par;, fi:R/Z/.,OA,ALPHAA.r,rTAA,r>AfWA.VOA,  ZB, QB.ALPHAe,  BETAS, 

1 GAW1AB , VOB , X AL  PHA  , REL  AX . F 1 , E2  , A 
COMMON/SFACTR/COSGTAU( 181 ) ,S!NGTAU.  181 ) 

C0riM0N/TABLE/IT(3,181  ),NT(11 ,11 ,11 ) 

C0Mi‘-10N/0IFF7AR(  181 ) ,AI  (181 ) 

COMMON/CHGDEN/Cn( 181 ) ,CI ( 181 ) 

C0KM0N/PRMTR2/ZC,QC,ALPHAC,BETAC,GAMMAC,V0C,2D,qD,AUPHA0,B£TAD, 

IGAMMAD.VOD.Xl ,Y1 
C 
C 
C 

AR(1)=0. 

AI(1)=0. 

C 

C=2. 5465/A 
C 

C1A=(1.-X1)*(ZA-QA)+X1*(ZC'QC) 

C1B=(1.-Y1)*(ZB-QB)+Y1*(ZD-Q0) 

C 

C2A=-19.7392*((1.-X1)/(BETAA*A)**2+X1/(BETAC*A)**2) 

C2B=-19.7392*((1.-Y1)/(BETAC*A)**2+Y1/(6ETAD*A)*»2) 

C 

C 

DO  200  1=2,181 

6SQ=1T(1,I)**2+IT(2,I)**2+IT(3  I)**2 
C 

RH0A=C1A*EXP(C2A*GSQ) 

RilOB=ClB*£XP(C2B*GSQ) 

C 

AR(  I ) =-  { C/GSQ) * ( CR ( I ) - (RHOB+P.HOA ) *COSGTAU(  I ) ) 
AI(I)=-(C/GSQ)*(CI(I)-(RHOB-RHOA)*SINGTAU(I)) 

C 

C 

200  CONTINUE 
C 
C 

RETURN 

END 
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SUBROUTINE  CLUSTER  (A,XA,RMAX,X1 ,Y1 ,Z) 

C 

C 

C THIS  SUBROUTINE  EVALUATES  A CLUSTER  SUMMATION  OVER  A RADIUS 

C RMAX  ABOUT  A POINT  X OF  FUNCTIONS  FTN  LOCATED  ON  ZINC-BLENDE 

C LATTICE  POINTS. 

C A IS  THE  CUBIC  LATTICE  CONSTANT. 

C 

C 

DIMENSION  X(3),XA(3),Y{3),TAU(3) 

C 

C 

LOGICAL  IX. lY 
C 

1X=X1.NE.0. 

IY=Yl.NE.O. 

C 

C 

Z=0. 

C 

DO  50  1 = 1 .3 
X(1)=A*XA{I) 

50  TAU0)  = .125*A 
C 

RMAX2=RMAX**2 

C=A/2. 

f'lAX=8 

MIN=-8 

C 

C 

C CONSTRUCT  THE  LATTICE  VECTORS. 

C 

DO  300  I1=MIN,MAX 
DO  300  I2=MIN,MAX 
DO  300  I3=MIN,MAX 
C 

Y(1)=C*(I2+I3) 

Y(2  =C*(n  + I3) 

Y(3)=C*{IUI2) 

C 

C A ATOM 

C 

R2=0. 

DO  100  1=1 ,3 

100  R2=R2+(Y(I)+TAU(I)-X(I))**2 
IF  (R2-KMAX2)  120,120,150 
120  R=SQRT(R2) 

Z=Z+(1,-X1)*VAT0MA(R) 

C 

IF  (IX)  Z=Z+X1*VAT0MC(R) 

C 

C B ATOM 

C 

150  R2=0. 

DO  170  1=1,3 

170  R2=R2+(Y(I)-TAU(I)-X(1))**2 
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IF  (R2-1!MAX2)  200,200,300 
200  R=SQRT(R2) 

Z=Z+(1.-Y1)*VAT0MB(R) 

IF(IY)  Z=Z+Yl*VAT0Hn(R) 


300  CONTINUE 
C 
C 

RETURN 

END 


L 


. t, 


i 
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FUNCTION  ERFC  (X) 

THIS  FUNCTION  EVALUATES  THE  COMPLEMENTARY  ERROR  FUNCTION. 


T=1./(1.+.47047*X) 

ERFC=T*  (,  3480242U*  {-.  0958798+T* , 7478556 ) )*EXP  ( -X**2 ) 

RETURN 

END 


FUNCTION  VATOrW  (R) 


C0MH0N/PARAriTR/?A,QA,ALPtiAA,BETAA,r,AMHAA,V0A,2B,QB,ALPHAR,B£:TAB. 
IGAMMAB.VOB.XALPHA.RELAX.El ,E2,A 

X=.707107*BETAA*R 

VAT0MA=(2./R)*(QA*EXP(-ALPHAA*R)+(ZA-QA)*ERFC(X)) 

RETURN 

END 


FUNCTION  VATO.MB  (R) 


COMMON/PARAMTR/ZA,OA,ALPHAA,BETAA,GAMMAA.VOA,ZC,QB.ALPHAB,BETAB, 
IGAMMAB.VOB.XALPHA.RELAX.El ,E2,A 

X=.70710Z*BETAB*R 

VAT0MB=(2./R)*(QB*EXP(-ALPHAB*R)+(ZB-QB)*ERFC(X)) 

RETURN 

END 


FUNCTION  VATOnC  (R) 


COMMON/PRf1TR2/ZC,QC,ALPHAC,BETAC.GAririAC,VOC,ZD,QD.ALPHAD,BETAD. 
IGAMMAD.VOD.Xl ,Y1 

X=.707107*BETAC*R 

VAT0.MC=(2./R)*(QC*EXP(-ALPHAC*R)  + (ZC-QC)*ERFC(X)) 

RETURN 

END 


FUNCTION  VATOriD  (R) 


COMMON/PRMTR2/ZC ,QC , ALPHAC , BETAC , GAMMAC , VOC , ZD ,QD , ALPHAO, BETAD , 
IGAMMAD.VOD.Xl  .Y1 

X=.707107*BETAD*R 

VAT0MD=(2./R)*(QD*EXP(-ALPHAD*R)+(ZD-QD)*ERFC(X)) 

RETURN 

END 


SUBROUTINE  FSPT  (N.X.AR.AI ,Z) 


C0MM0N/TABLE/1T(3.181).NT( 11,11 ,11) 

DIMENSION  AP(fi),AI(N) 

DIMENSION  X(3) 


DO  200  1 = 1 ,N 
GX=0. 

DO  100  J=l,3 
100  GX=r.X+lT(J,I)*X{J) 

GX=6.28319*GX 

2=Z+AR(I)*C0S(GX)-AI(I)*SIN(GX) 
200  CONTINUE 


f~)  C">  o o o oooooo  ooo 


SUBROUTIfit  FR5BM  (N,XA,XB,AP,AI  ,Z) 

C 

C 

C 

COMMON/ TABLE/ 1 T ( 3 , 1 81 ) .NT ( 11  J 1 . 1 1 ) 

C 

C 

DIMENSION  XA(3),XB(3),DELX(3) 

DIMENSION  AR(181),A1(IG1) 

DIMENSION  ?{N) 

C 

C 

C=6. 28319 

CALCULATE  DELX 
DO  50  1-1,3 

50  DELX(I)-(XB(I)-XA(I))/N 
DO  FOR  ALL  G 
DO  400  J=1 .229 

INITIALIZE  COSGX.SINGX.COSGDX.SINGDX 

GX=0. 

GDX-0. 

00  100  1-1,3 
GX-GX+IT(I,J)*XA(I) 

100  GDX=GDX+n(I,J)*DELX(I) 

C0SGX=C0S(C*GX) 

SINGX  = SIN(C*r.X) 

COSGDX-COS(f’GDX) 

S1NGDX=SIN(C*GDX) 

DO  FOR  ALL  X 

DO  300  K=1  ,N 
C 

COST-COSGX 

SINT-SINGX 

C 

COSGX=COST*COSGDX-SINT*SINGOX 

S1NGX=C0ST*SINGDX+SINT*C0SGDX 

C 

Z(K)=Z(K)+AR(J)*C0SGX-AI(J)*S1NGX 

C 

C 

300  CONTINUE 
C 

400  CONTINUE 
C 

RETURN 

END 


c 

c 

c 


SUBROUTINE  PLOTR  (Y.Z£RO,R/>:iC£) 


DIMENSION  Y(100),.NAME(8) 

DIMENSION  ICHAR(IOO) 

C 

C 

READ  100,(NAM£(I),I=1 .8) 

100  FORMAT  (8A10) 

C 

PRINT  110,(NAM£(I).!-1 .8) 
no  FORMAT  (*1*,10X,SA10,///) 

C 

c 

D£l.Y=RANr,£/50. 

IFLG-50 

C 

C 

c 

DO  300  INDX=1 .51 
C 

M=51-INDX 

YM=M»DELY+ZERO 

YMIN=Yf1-.5*DELY 

YMAX=YM+.5*DELY 

C 

C 

DO  220  K=l,100 
C 

ICHAR{K)=0000055 

C 

IF  ((Y(K).LE.YMAX).AND.(Y(K).GT.YMIN))  ICHAR( K)=0000047 
C 

220  CONTINUE 
C 
C 

IF  (M.EQ.IFLG)  GO  TO  230 
C 

PRINT  225,(ICHAR(K),K=1,100) 

225  FORMAT  (15X,*I*,100R1 ) 

C 

GO  TO  300 
C 

c 

230  PRINT  235,YM,{ICHAR(K),K=1 ,100) 

235  FORMAT  (5X,F10.2,*I*.100R1 ) 

IFLG=IFLG-5 

C 

C 

300  CONTINUE 
C 
C 
C 

PRINT  310 

310  FORMAT  (15X,20(*I ....*) .*!•) 


